-
Notifications
You must be signed in to change notification settings - Fork 15
/
stsubs.F90
1292 lines (1252 loc) · 39 KB
/
stsubs.F90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!Standard procedures often used by me.
!AUTHOR: Dmitry I. Lyakh (Liakh): [email protected]
!REVISON: 2020/05/29
!Copyright (C) 2005-2022 Dmitry I. Lyakh (Liakh)
!LICENSE: BSD 3-Clause
module stsubs
use, intrinsic:: ISO_C_BINDING
implicit none
private
!Parameters:
logical, private:: VERBOSE=.false. !verbosity for errors
real(8), parameter, public:: PI=3.14159265358979d0 !PI constant
real(8), parameter, public:: BOHR=0.529177249d0 !Bohrs to Angstroms conversion factor
!Procedures:
public:: alphanumeric!checks if the character is alphanumeric
public:: alphanumeric_underscore!checks if the character is alphanumeric or underscore
public:: alphanumeric_string !checks if the string is alphanumeric (with underscores allowed)
public:: array2string!converts a character*1 array into a string
public:: byte_chksum !returns the raw byte check sum for a given memory segment
public:: cap_ascii !makes all English letters capital
public:: charnum !converts a number given as a string to real*8 and integer numbers
public:: clone_object!clones a Fortran object
public:: crash !causes a crash
public:: create_line !creates a table line in .txt format with ; separator
public:: dumb_work !performs a dumb work on one or two arrays producing a third one
public:: dump_bytes !dumps byte values for a given memory segment
public:: icharnum !converts a number given as a string to the integer number
public:: ifcl !calculates factorial
public:: is_it_letter!checks if the character is an ASCII letter
public:: is_it_number!checks if the character is an ASCII number
public:: is_nan_c4 !checks whether a complex(4) number is NaN
public:: is_nan_c8 !checks whether a complex(8) number is NaN
public:: is_nan_r4 !checks whether a real(4) number is NaN
public:: is_nan_r8 !checks whether a real(8) number is NaN
public:: itrsign !determines a sign of a given transposition
public:: longnumchar !converts a long integer number to the character representation
public:: markchf !counts how many non-blank fields a string contains
public:: matmat !multiplies matrix on matrix
public:: matmatt !multiplies matrix on transposed matrix
public:: mattmat !multiplies transposed matrix on matrix
public:: mattmatt !multiplies transposed matrix on transposed matrix
public:: name_hash !returns a hash-mask for a given string
public:: nocomment !removes comments from a line (!,#)
public:: normv !normalizes a vector
public:: nospace !removes blanks from a string
public:: not_a_number!checks whether a given string solely contains a decimal number
public:: numchar !converts an integer number to character representation
public:: printl !prints a line of characters
public:: rand_bool !returns a random logical
public:: rand_int4 !returns a random integer4 from the given range
public:: rand_real8 !returns a random real8
public:: rand_str !returns a random string
public:: rots !rotates an array of points in a 3d-space
public:: size_of !scalar data type size in bytes (Fortran 2008)
public:: small_ascii !makes all English letters small
public:: str_cmp !compares two strings
public:: string2array!converts a string into a character*1 array
public:: strsearch !searches a given fragment in a string.
public:: symbol_type !function which tells you whether the given symbol is a space/tab (0), number (1), letter (2), or other (-1)
public:: tpause !pause for a given number of seconds
public:: valchar !converts a real number to the string of characters
public:: wait_delay !pause for a given time
public:: wait_press !subroutine makes pause until user presses a key
public:: wr_mat_in !writes a matrix of integers to the screen
public:: wr_mat_in8 !writes a matrix of integer8's to the screen
public:: wr_mat_sp !writes a matrix of single precision elements to the screen
public:: wr_mat_dp !writes a matrix of double precision elements to the screen
public:: wr_mat_dc !writes a matrix of double complex elements to the screen
public:: wr_vec_sp !writes a vector of single precision elements to the screen
public:: wr_vec_dp !writes a vector of double precision elements to the screen
contains
!----------------------------------------
logical function alphanumeric(ch)
!Returns TRUE if the character is ASCII alphanumeric, FALSE otherwise.
!A space is not considered alphanumeric.
implicit none
character(1), intent(in):: ch
alphanumeric=.FALSE.
if(is_it_letter(ch).gt.0) then
alphanumeric=.TRUE.
else
if(is_it_number(ch).ge.0) alphanumeric=.TRUE.
endif
return
end function alphanumeric
!---------------------------------------------------
logical function alphanumeric_underscore(ch)
!Returns TRUE if the character is ASCII alphanumeric or underscore, FALSE otherwise.
!A space is not considered alphanumeric.
implicit none
character(1), intent(in):: ch
alphanumeric_underscore=.FALSE.
if(is_it_letter(ch).gt.0) then
alphanumeric_underscore=.TRUE.
else
if(is_it_number(ch).ge.0) then
alphanumeric_underscore=.TRUE.
else
if(ch.eq.'_') alphanumeric_underscore=.TRUE.
endif
endif
return
end function alphanumeric_underscore
!------------------------------------------------
logical function alphanumeric_string(str)
!Returns TRUE if the string only contains ASCII alphanumeric + underscore,
!FALSE otherwise. Spaces are not considered alphanumeric_.
!An empty string is not considered alphanumeric_ either.
implicit none
character(*), intent(in):: str
integer:: i,l
alphanumeric_string=.TRUE.
l=len(str)
if(l.gt.0) then
do i=1,l
if(.not.alphanumeric_underscore(str(i:i))) then
alphanumeric_string=.FALSE.
exit
endif
enddo
else
alphanumeric_string=.FALSE.
endif
return
end function alphanumeric_string
!------------------------------------------------
subroutine array2string(str,ar1,arl,ierr)
!Converts a CHARACTER(1) array into a Fortran string.
implicit none
character(*), intent(out):: str
integer, intent(in):: arl
character(1), intent(in):: ar1(1:arl)
integer i,ierr
ierr=0
if(arl.le.len(str)) then
do i=1,arl; str(i:i)=ar1(i); enddo
else
ierr=1
endif
return
end subroutine array2string
!----------------------------------------------------
function byte_chksum(cptr,csize) result(csum)
!Returns the raw byte check sum for a given memory segment.
implicit none
integer(8):: csum !out: check sum
type(C_PTR), intent(in):: cptr !in: base C pointer
integer, intent(in):: csize !in: size in bytes
integer(1), pointer:: iptr(:)
integer:: i
csum=0_8
call c_f_pointer(cptr,iptr,shape=(/csize/))
do i=1,csize
if(iptr(i).ge.0) then
csum=csum+int(iptr(i),8)
else
csum=csum+(256_8+int(iptr(i),8))
endif
enddo
return
end function byte_chksum
!--------------------------------
subroutine cap_ascii(str)
!Capitalizes all small English letters in string <str>
implicit none
character(*), intent(inout):: str
integer i,k,l,ia1,ia2,ia3
ia1=iachar('a')
ia2=iachar('z')
ia3=iachar('A')
l=len_trim(str)
do k=1,l
i=iachar(str(k:k))
if(i.ge.ia1.and.i.le.ia2) str(k:k)=achar(ia3+(i-ia1))
enddo
return
end subroutine cap_ascii
!---------------------------------
subroutine charnum(A,DN,N)
!Converts a number given as a string (<A>) to real and integer forms.
!If <A> is not a number, then it will set {int(DN)=-2, N=0} (not a number tag).
!A real number in <A> can use both "." and "," for the decimal point designation.
!If <A> contains a number at the beginning separated by a blank/tab from the
!following part of the string, only this number will be converted while
!the rest of the string will be ignored (if ALLOW_END_GARBAGE set to TRUE).
!Note that <N> is only guaranteed to be equal to the original (string) number
!when that number is an integer written either as an integer or as a scientific
!number with a bare (D|d|E|e) or with a trvial (D0|d0|E0|e0). In all other cases,
!<N> is not required (or cannot) represent the original (string) number.
implicit none
character(*), intent(in):: A !in: string containing a symbolic number (integer or real)
integer, intent(out):: N !out: integer number (rounded off for reals)
real(8), intent(out):: DN !out: real number
logical, parameter:: ALLOW_END_GARBAGE=.true.
real(8):: DSN,DP,DPS
integer:: I,J,K,K1,K2,L,M,KF,LA,IA0,IA9,ISN,ADP,ierr
ierr=0; N=0; DN=0d0
LA=len_trim(A)
if(LA.gt.0) then
!Check the sign of the number:
IA0=iachar('0'); IA9=iachar('9')
DSN=+1d0; ISN=+1 !default sign is positive
DP=0d0; DPS=+1d0 !power (for scientific D and E formats)
M=1; KF=0; ADP=0
K1=1; do while(K1.lt.LA.and.(A(K1:K1).eq.' '.or.iachar(A(K1:K1)).eq.9)); K1=K1+1; enddo !skip initial blanks
if(A(K1:K1).eq.'+') then !positive number
if(K1.lt.LA) then; K1=K1+1; else; ierr=1; endif
elseif(A(K1:K1).eq.'-') then !negative number
if(K1.lt.LA) then; DSN=-1d0; ISN=-1; K1=K1+1; else; ierr=2; endif
endif
else
ierr=3
endif
!Process the number itself:
if(ierr.eq.0) then
cloop: do K=K1,LA
I=iachar(A(K:K))
if(IA0.le.I.and.I.le.IA9) then !decimal digit
KF=max(KF,1) !decimal digit found
if(M.eq.1) then !integer part
N=N*10+(I-IA0)
DN=DN*10d0+real(I-IA0,8)
else !fractional part
ADP=ADP-1
DN=DN+real(I-IA0,8)*(10d0**real(ADP,8))
if(I.gt.IA0) KF=2 !non-zero fractional part
endif
elseif(A(K:K).EQ.'.'.or.A(K:K).EQ.',') then !decimal dot/comma detected
if(M.eq.1) then; M=M+1; else; ierr=4; exit cloop; endif
elseif(A(K:K).EQ.'D'.or.A(K:K).EQ.'d'.or.A(K:K).EQ.'E'.or.A(K:K).EQ.'e') then !power symbol
if(KF.ne.0) then
L=K+1
if(L.gt.LA) exit cloop !bare D means D0 by default (the same for d,E,e)
J=iachar(A(L:L))
if(A(L:L).eq.'+') then
if(L.lt.LA) then; L=L+1; else; ierr=5; exit cloop; endif
elseif(A(L:L).eq.'-') then
if(L.lt.LA) then; DPS=-1d0; L=L+1; else; ierr=6; exit cloop; endif
elseif(J.lt.IA0.or.J.gt.IA9) then
ierr=7; exit cloop !not a number case
endif
do K2=L,LA
J=iachar(A(K2:K2))
if(J.ge.IA0.and.J.le.IA9) then
if(J.gt.IA0) KF=2 !non-zero power value
DP=DP*10d0+real(J-IA0,8)
else
ierr=8; exit cloop !not a number case
endif
enddo
DP=DPS*DP
exit cloop
else
ierr=9; exit cloop
endif
elseif(A(K:K).eq.' '.or.I.eq.9) then !blank/tab
if(KF.eq.0.or.(.not.ALLOW_END_GARBAGE)) ierr=10
exit cloop
else !unexpected character
ierr=11; exit cloop !not a number case
endif
enddo cloop !next character A(K:K)
endif
!Result:
if(ierr.eq.0.and.KF.ne.0) then !number succesfully translated
DN=DSN*(DN*(10d0**DP))
if(KF.gt.1) then
N=int(DN)
else
N=ISN*N
endif
else !not a number case
DN=-2d0; N=0 !{DN=-2d0 & N=0}: not a number case
endif
return
end subroutine charnum
!-------------------------------------------------------
function clone_object(object,ierr) result(clone)
!Clones a Fortran object: Deep copy for allocatable components,
!pointer association for pointer components.
implicit none
class(*), pointer:: clone !out: clone
class(*), intent(in), target:: object !in: object
integer, intent(out), optional:: ierr !out: error code
character(256):: errmesg
integer:: errc
clone=>NULL(); errc=0
allocate(clone,SOURCE=object,STAT=errc,ERRMSG=errmesg)
if(errc.ne.0) then
write(*,*)'#ERROR(stsubs:clone_object): Fortran sourced allocate() failed: '//errmesg
if(errmesg(1:39).eq.'Attempt to allocate an allocated object') then !debug
write(*,*)'If you see status = F below, you are likely experiencing a gfortran bug! Switch to GCC/8.0.0+'
write(*,*)'Object (pointer) association status = ',associated(clone)
endif
clone=>NULL()
endif
if(present(ierr)) ierr=errc
return
end function clone_object
!-------------------------
subroutine crash()
integer:: ax,dx
write(*,'("Initiating a crash by dividing by zero ...")')
ax=1; dx=0; ax=ax/dx; write(*,*) ax
return
end subroutine crash
!------------------------------------------------------------------
subroutine create_line(talen,tabln,mantl,ctrls,numar,symar)
!The subroutine creates a table line in .txt format with some separator.
!INPUT:
! 1) mantl - number of mantissa digits;
! 2) ctrls - control sequence of commands;
! 3) symar - symbolic array of data (pool of symbolic entries);
! 4) numar - real(8) array of data (pool of numeric entries);
!OUTPT:
! 1) tabln - table line (string of characters);
! 2) talen - length of <tabln>.
implicit none
!-----------------------------------------------------------
character(1), parameter:: sep=';' !default separator
!-----------------------------------------------------------
integer i,j,k,l,m,n,k1,k2,k3,k4,kf
integer, intent(in):: mantl
integer, intent(out):: talen
character(*), intent(out):: tabln
integer, intent(in):: ctrls(*)
real(8), intent(in):: numar(*)
character(*), intent(in):: symar(*)
character(256):: os
talen=0
l=1 !instruction pointer
do
select case(ctrls(l))
case(0) !end of control sequence
exit
case(1) !insert a number of separators
do k=1,ctrls(l+1)
talen=talen+1
tabln(talen:talen)=sep
enddo
l=l+2
case(2) !insert symbolic entries separated by separators
do k=ctrls(l+1),ctrls(l+2)
m=len_trim(symar(k))
tabln(talen+1:talen+m+1)=symar(k)(1:m)//sep
talen=talen+m+1
enddo
l=l+3
case(3) !insert numeric entries separated by separators
do k=ctrls(l+1),ctrls(l+2)
call valchar(numar(k),mantl,m,os)
tabln(talen+1:talen+m+1)=os(1:m)//sep
talen=talen+m+1
enddo
l=l+3
case default
write(*,*)'#ERROR(STSUBS::CREATE_LINE): invalid control command: ',ctrls(l),l
stop
end select
enddo
return
end subroutine create_line
!-------------------------------------------
subroutine dumb_work(arr0,arr1,arr2)
!Performs a dumb work on one or two arrays, producing a third one.
implicit none
real(8), intent(inout), contiguous:: arr0(0:)
real(8), intent(in), contiguous:: arr1(0:)
real(8), intent(in), contiguous, optional:: arr2(0:)
integer:: vol0,vol1,vol2,i,j,k
real(8):: val,ax
vol0=size(arr0); vol1=size(arr1)
if(vol0.le.0.or.vol1.le.0) return
if(present(arr2)) then
vol2=size(arr2)
!$OMP PARALLEL DO DEFAULT(SHARED) SCHEDULE(GUIDED) PRIVATE(i,j,k,val,ax)
do i=0,vol0-1
ax=0d0
do j=mod(i,vol1),vol1-1
val=arr1(j)
do k=mod(i,vol2),vol2-1
ax=ax+arr2(k)*val
enddo
enddo
arr0(i)=ax
enddo
!$OMP END PARALLEL DO
else
!$OMP PARALLEL DO DEFAULT(SHARED) SCHEDULE(GUIDED) PRIVATE(i,j)
do i=0,vol0-1
arr0(i)=0d0
do j=mod(i,vol1),vol1-1
arr0(i)=arr0(i)+arr1(j)*arr1(j)
enddo
enddo
!$OMP END PARALLEL DO
endif
return
end subroutine dumb_work
!--------------------------------------------------
subroutine dump_bytes(cptr,csize,file_name)
!Dumps byte values for a given memory segment.
implicit none
type(C_PTR), intent(in):: cptr !in: base C pointer
integer, intent(in):: csize !in: size in bytes
character(*), intent(in), optional:: file_name !in: output file name (defaults to screen)
integer(1), pointer:: iptr(:)
integer:: i,devo
if(present(file_name)) then
devo=666
open(devo,file=file_name,FORM='FORMATTED',STATUS='UNKNOWN')
else
devo=6
endif
write(devo,'()'); write(devo,*) '### Memory dump for address ' !,cptr
call c_f_pointer(cptr,iptr,shape=(/csize/))
do i=1,csize
if(iptr(i).ge.0) then
write(devo,'("Offset ",i10,": ",i4)') i-1,int(iptr(i),4)
else
write(devo,'("Offset ",i10,": ",i4)') i-1,(256+int(iptr(i),4))
endif
enddo
if(devo.ne.6) close(devo)
return
end subroutine dump_bytes
!--------------------------------------
integer function icharnum(L,OS)
!Converts an integer number OS(1:L) given as a string into INTEGER.
!Resets <L> to 0 in case of error.
implicit none
integer, intent(inout):: L !inout: length of the string (resets to ZERO in case of error!)
character(*), intent(in):: OS !in: OS(1:L): positive or negative number as a string
integer:: K,M,N,IA0
IA0=iachar('0')
if(L.gt.0) then
if(OS(1:1).eq.'-') then !negative
M=-2
elseif(OS(1:1).eq.'+') then !explicitly positive
M=2
else !implicitly positive
M=1
endif
icharnum=0
do K=abs(M),L
N=iachar(OS(K:K))-IA0
if(N.ge.0.and.N.le.9) then
icharnum=icharnum*10 + N
else
if(VERBOSE) write(*,*)'#ERROR(STSUBS::ICHARNUM): invalid character number: '//OS(1:L)
L=0 !error
return
endif
enddo
if(M.eq.-2) icharnum=-icharnum
else
if(VERBOSE) write(*,*)'#ERROR(STSUBS::ICHARNUM): string of non-positive length: ',L
L=0 !error
endif
return
end function icharnum
!-------------------------------
integer function ifcl(N)
!Returns the factorial of a non-negative N.
implicit none
integer N,K
ifcl=1
do K=2,N
ifcl=ifcl*K
enddo
return
end function ifcl
!----------------------------------------
integer function is_it_letter(ch)
!If character <ch> is a letter, returns 1 for lower-case, 2 for upper-case,
!otherwise returns 0.
implicit none
character(1), intent(in):: ch
integer:: i
i=iachar(ch)
if(i.ge.iachar('a').and.i.le.iachar('z')) then
is_it_letter=1
elseif(i.ge.iachar('A').and.i.le.iachar('Z')) then
is_it_letter=2
else
is_it_letter=0
endif
return
end function is_it_letter
!----------------------------------------
integer function is_it_number(ch)
!If character <ch> is a number, returns its value as an integer,
!otherwise a negative integer is returned.
implicit none
character(1), intent(in):: ch
integer:: i
i=iachar(ch)
if(i.ge.iachar('0').and.i.le.iachar('9')) then
is_it_number=i-iachar('0')
else
is_it_number=-1
endif
return
end function is_it_number
!--------------------------------------
logical function is_nan_c4(fpn)
!Returns TRUE if the <fpn> number is NaN.
complex(4), intent(in):: fpn
is_nan_c4 = (fpn.ne.fpn)
return
end function is_nan_c4
!--------------------------------------
logical function is_nan_c8(fpn)
!Returns TRUE if the <fpn> number is NaN.
complex(8), intent(in):: fpn
is_nan_c8 = (fpn.ne.fpn)
return
end function is_nan_c8
!--------------------------------------
logical function is_nan_r4(fpn)
!Returns TRUE if the <fpn> number is NaN.
real(4), intent(in):: fpn
is_nan_r4 = (fpn.ne.fpn)
return
end function is_nan_r4
!--------------------------------------
logical function is_nan_r8(fpn)
!Returns TRUE if the <fpn> number is NaN.
real(8), intent(in):: fpn
is_nan_r8 = (fpn.ne.fpn)
return
end function is_nan_r8
!--------------------------------
subroutine itrsign(N,ITR)
!Reorders a given permutation into an ascending order and returns the permutation sign in ITR(0).
!Bubble algorithm (only for small permutations).
implicit none
integer, intent(inout):: ITR(0:*) !inout: ITR(1:N) is the permutation, ITR(0) is its initial/final sign
integer, intent(in):: N !in: length of the permutation
integer:: K,L
K=1
do while(K.lt.N)
if(ITR(K).gt.ITR(K+1)) then
L=ITR(K); ITR(K)=ITR(K+1); ITR(K+1)=L
ITR(0)=-ITR(0)
if(K.gt.1) then
K=K-1
else
K=2
endif
else
K=K+1
endif
enddo
return !RETURNS an ascending sequence of indices with a sign in ITR(0)
end subroutine itrsign
!----------------------------------------
subroutine longnumchar(I,IOSL,OS)
!Converts a long integer into a string.
implicit none
integer(8), intent(in):: I !in: long integer to be converted into a string
integer, intent(out):: IOSL !out: string length
character(*), intent(inout):: OS !out: OS(1:IOSL) is the string
integer(8) K
integer L,K1,K2
character(1) A(0:9),CH
data A/'0','1','2','3','4','5','6','7','8','9'/
if(I.lt.0) then
OS(1:1)='-'; K=-I; IOSL=1
elseif(I.eq.0) then
OS(1:1)='0'; IOSL=1
return
else
K=I; IOSL=0
endif
L=IOSL
do while(K.NE.0)
K1=mod(K,10_8); K=K/10_8
L=L+1; OS(L:L)=A(K1)
enddo
K1=L-IOSL; if(mod(K1,2).eq.1) K1=K1-1; K1=K1/2
do K2=1,K1
CH=OS(IOSL+K2:IOSL+K2)
OS(IOSL+K2:IOSL+K2)=OS(L+1-K2:L+1-K2)
OS(L+1-K2:L+1-K2)=CH
enddo
IOSL=L
return !OS(1:IOSL) - the number as a string (the sign included if the number is negative)
end subroutine longnumchar
!-----------------------------------
subroutine markchf(STR,N,MF)
!Counts the number of fields in the string <STR> (separated by spaces) and returns
!the positions of all fields found (a field is a sequence of characters without blanks).
implicit none
character(*), intent(in):: STR !in: string
integer, intent(out):: N !out: N is the number of fields in the string
integer, intent(out):: MF(2,1:*) !out: MF(1:2,X) is the beginning and end positions of the field #X(1<=X<=N)
integer:: K,LSTR
LSTR=len_trim(STR)
N=0; K=1
do while(K.le.LSTR)
if(STR(K:K).ne.' '.and.iachar(STR(K:K)).ne.9) then !skip TABS also (ASCII code 9)
N=N+1; MF(1:2,N)=(/K,0/)
do while(K.le.LSTR)
if(STR(K:K).eq.' '.or.iachar(STR(K:K)).eq.9) exit
K=K+1
enddo
MF(2,N)=K-1
else
K=K+1
endif
enddo
return
end subroutine markchf
!---------------------------------------
subroutine matmat(sa,sb,m,a,b,c)
!Matrix-matrix multiplication: C(sa,sb)=A(sa,m)*B(m,sb)+C(sa,sb)
implicit none
integer, intent(in):: sa,sb,m
real(8), intent(in):: a(1:sa,1:m),b(1:m,1:sb)
real(8), intent(inout):: c(1:sa,1:sb)
integer j,k
do j=1,sb
do k=1,m
c(1:sa,j)=c(1:sa,j)+a(1:sa,k)*b(k,j)
enddo
enddo
return
end subroutine matmat
!----------------------------------------
subroutine matmatt(sa,sb,m,a,b,c)
!Matrix-matrix multiplication (transposed 2nd argument): C(sa,sb)=A(sa,m)*B(sb,m)+C(sa,sb)
implicit none
integer, intent(in):: sa,sb,m
real(8), intent(in):: a(1:sa,1:m),b(1:sb,1:m)
real(8), intent(inout):: c(1:sa,1:sb)
integer j,k
do j=1,sb
do k=1,m
c(1:sa,j)=c(1:sa,j)+a(1:sa,k)*b(j,k)
enddo
enddo
return
end subroutine matmatt
!----------------------------------------
subroutine mattmat(sa,sb,m,a,b,c)
!Matrix-matrix multiplication (transposed 1st argument): C(sa,sb)=A(m,sa)*B(m,sb)+C(sa,sb)
implicit none
integer, intent(in):: sa,sb,m
real(8), intent(in):: a(1:m,1:sa),b(1:m,1:sb)
real(8), intent(inout):: c(1:sa,1:sb)
integer i,j,k
real(8) val
do j=1,sb
do k=1,sa
val=c(k,j)
do i=1,m; val=val+a(i,k)*b(i,j); enddo !unroll this loop
c(k,j)=val
enddo
enddo
return
end subroutine mattmat
!-----------------------------------------
subroutine mattmatt(sa,sb,m,a,b,c)
!Matrix-matrix multiplication (both arguments are transposed): C(sa,sb)=A(m,sa)*B(sb,m)+C(sa,sb)
implicit none
integer, intent(in):: sa,sb,m
real(8), intent(in):: a(1:m,1:sa),b(1:sb,1:m)
real(8), intent(inout):: c(1:sa,1:sb)
integer j,k
do j=1,sb
do k=1,m
c(1:sa,j)=c(1:sa,j)+a(k,1:sa)*b(j,k)
enddo
enddo
return
end subroutine mattmatt
!--------------------------------------------------
integer function name_hash(max_hash,sl,str)
!This function returns a hash mask (0..max_hash) for a given string.
!INPUT:
! - max_hash - limit of the hash range (parameter);
! - str(1:sl) - string;
!OUTPUT:
! - name_hash - hash mask [0..max_hash].
implicit none
integer, intent(in):: max_hash,sl
character(*), intent(in):: str
integer:: i
name_hash=0
do i=1,min(sl,len(str))
name_hash=name_hash+iachar(str(i:i))
enddo
name_hash=mod(name_hash,max_hash+1)
return
end function name_hash
!----------------------------------
subroutine nocomment(l,str)
!This subroutine removes comments (!,#) from the line given in str(1:l).
!After that l is a modified length.
implicit none
integer, intent(inout):: l
character(*), intent(in):: str
integer:: i
do i=1,min(l,len(str))
if(str(i:i).eq.'!'.or.str(i:i).eq.'#') then
l=i-1
exit
endif
enddo
return
end subroutine nocomment
!-------------------------------
subroutine normv(n,vn,v)
!Substitutes the original vector <v> with its normalized form.
implicit none
integer, intent(in):: n !in: dimension of the vector
real(8), intent(out):: vn !out: vector length
real(8), intent(inout):: v(1:n) !inout: initial vector will be replaced by the normalized one
integer:: i
real(8):: val
vn=0d0
do i=1,n
vn=vn+v(i)*v(i)
enddo
if(vn.ne.0d0) then
vn=dsqrt(vn); val=1d0/vn
v(1:n)=v(1:n)*val
endif
return
end subroutine normv
!------------------------------
subroutine nospace(A,L)
!Removes spaces/tabs from a string.
implicit none
character(*), intent(inout):: A !inout: string --> reduced string
integer, intent(out):: L !out: the length of the string with all spaces/tabs removed
integer:: J,M
M=len_trim(A); L=0
do J=1,M
if(A(J:J).ne.' '.and.iachar(A(J:J)).ne.9) then
L=L+1; A(L:L)=A(J:J)
endif
enddo
return
end subroutine nospace
!-----------------------------------------
logical function not_a_number(str)
!Returns TRUE if string <str> is a number, false otherwise.
!Signed integer, decimal, and scientific decimal formats are supported.
!The decimal point is ".". For scientific decimals, the exponent can be either of {D|d|E|e}.
implicit none
character(*), intent(in):: str
integer:: dot_present,exp_present,i,l,n
l=len_trim(str)
if(l.gt.0) then
not_a_number=.false.; dot_present=0; exp_present=0; n=0
if(str(1:1).eq.'-'.or.str(1:1).eq.'+') then; i=2; else; i=1; endif
do while(i.le.l)
if(iachar(str(i:i)).lt.iachar('0').or.iachar(str(i:i)).gt.iachar('9')) then
if(str(i:i).eq.'.') then !decimal dot
if(dot_present.eq.0.and.exp_present.eq.0) then; dot_present=1; else; not_a_number=.true.; return; endif
elseif(str(i:i).eq.'D'.or.str(i:i).eq.'d'.or.str(i:i).eq.'E'.or.str(i:i).eq.'e') then
if(exp_present.eq.0.and.n.gt.0) then
exp_present=1
if(i.lt.l) then; if(str(i+1:i+1).eq.'-'.or.str(i+1:i+1).eq.'+') i=i+1; else; not_a_number=.true.; return; endif
else
not_a_number=.true.; return
endif
else
not_a_number=.true.; return
endif
else !decimal digit
n=n+1
endif
i=i+1
enddo
if(n.le.0) not_a_number=.true.
else !empty string is not a number
not_a_number=.true.
endif
return
end function not_a_number
!------------------------------------
subroutine numchar(I,IOSL,OS)
!Converts an integer into a string.
implicit none
integer, intent(in):: I !in: integer to be converted into a string
integer, intent(out):: IOSL !out: length of the string
character(*), intent(inout):: OS !out: OS(1:IOSL) is the string
integer K,K1,L
character(1) A(0:9),CH
data A/'0','1','2','3','4','5','6','7','8','9'/
if(I.lt.0) then
OS(1:1)='-'; K=-I; IOSL=1
elseif(I.eq.0) then
OS(1:1)='0'; IOSL=1
return
else
K=I; IOSL=0
endif
L=IOSL
do while(K.ne.0)
K1=mod(K,10); K=K/10
L=L+1; OS(L:L)=A(K1)
enddo
K1=L-IOSL; if(mod(K1,2).eq.1) K1=K1-1; K1=K1/2
do K=1,K1
CH=OS(IOSL+K:IOSL+K)
OS(IOSL+K:IOSL+K)=OS(L+1-K:L+1-K)
OS(L+1-K:L+1-K)=CH
enddo
IOSL=L
return !OS(1:IOSL) - the number as a string (the sign included if the number is negative)
end subroutine numchar
!------------------------------------
subroutine printl(FH,STR,ADV)
!Prints a string <str> as a single line into file #FH.
implicit none
integer, intent(in):: FH !in: file handle
character(*), intent(in):: STR !in: string to print
logical, intent(in), optional:: ADV !in: carriage return control (.false. - no carriage return)
character(32):: FRM
integer:: l,k
l=len_trim(STR)
if(l.gt.0) then
FRM(1:2)='(A'
call NUMCHAR(l,k,FRM(3:32))
k=3+k; FRM(k:k)=')'
if(present(ADV)) then
if(ADV) then
write(FH,FRM(1:k)) STR(1:l)
else
write(FH,FRM(1:k),advance='no') STR(1:l)
endif
else
write(FH,FRM(1:k)) STR(1:l)
endif
else
if(present(ADV)) then
if(ADV) write(FH,'()')
endif
endif
return
end subroutine printl
!----------------------------------------
function rand_bool() result(bool)
!Returns a random logical.
implicit none
logical:: bool !out: random logical
real(8):: rn
bool=.FALSE.
call random_number(rn)
if(rn.ge.5d-1) bool=.TRUE.
return
end function rand_bool
!---------------------------------------------------
function rand_int4(lower,upper) result(int4)
!Returns a random integer4 from the given range.
implicit none
integer(4):: int4 !out: random number
integer(4), intent(in):: lower !in: lower bound
integer(4), intent(in):: upper !in: upper bound
int4=nint(rand_real8(real(lower,8),real(upper,8)),4)
return
end function rand_int4
!-----------------------------------------------------
function rand_real8(lower,upper) result(real8)
!Returns a random real8 from the given range.
implicit none
real(8):: real8 !out: random number
real(8), intent(in):: lower !in: lower bound
real(8), intent(in):: upper !in: upper bound
call random_number(real8)
real8=lower+real8*(upper-lower)
return
end function rand_real8
!----------------------------------
subroutine rand_str(sl,str)
!Generates a random string.
implicit none
integer, intent(out):: sl !out: string length
character(*), intent(inout):: str !out: str(1:sl) is the random string
integer, parameter:: max_str_len=8192
integer, parameter:: buf_size=512
integer, parameter:: ascii_b=32,ascii_e=126,ascii_range=ascii_e-ascii_b+1
real(8) rnds(1:buf_size),rn
integer i,j,k,m,n
if(len(str).gt.0) then
k=min(len(str),max_str_len)
call random_number(rn); sl=int(rn*real(k,8))+1; if(sl.gt.k) sl=k
do i=0,sl-1,buf_size
n=min(sl-i,buf_size)
call random_number(rnds(1:n))
do j=1,n
m=int(rnds(j)*real(ascii_range,8))+ascii_b; if(m.gt.ascii_e) m=ascii_e
str(i+j:i+j)=achar(m)
enddo
enddo
else
sl=0
endif
return
end subroutine rand_str
!-----------------------------------------
subroutine rots(axis,alpha,n,pnts)
!This subroutine rotates <n> points in array <pnts>
!around the axis #axis in a 3d-space.
implicit none
integer i,j,k,l,m,n,k1,k2,k3,k4,ks,kf
integer axis
real(8) alpha,pnts(3,*),q(3),cosa,sina
sina=sin(alpha)
cosa=cos(alpha)
select case(axis)
case(1)
!X axis:
do k=1,n
q(1:3)=pnts(1:3,k)
pnts(1,k)=q(1)
pnts(2,k)=q(2)*cosa-q(3)*sina
pnts(3,k)=q(2)*sina+q(3)*cosa
enddo
case(2)
!Y axis:
do k=1,n
q(1:3)=pnts(1:3,k)
pnts(1,k)=q(1)*cosa-q(3)*sina
pnts(2,k)=q(2)
pnts(3,k)=q(1)*sina+q(3)*cosa
enddo
case(3)
!Z axis:
do k=1,n
q(1:3)=pnts(1:3,k)
pnts(1,k)=q(1)*cosa-q(2)*sina
pnts(2,k)=q(1)*sina+q(2)*cosa
pnts(3,k)=q(3)
enddo
case default
write(*,*)'#ERROR(STSUBS::ROTS): invalid axis number: ',axis
stop
end select
return
end subroutine rots
!-------------------------------------
integer function size_of(uarg) !Fortran 2008
!Returns size in bytes for an arbitrary scalar type (-1 means an error).
implicit none
class(*), intent(in), target:: uarg !in: scalar argument of any type/class
size_of=storage_size(uarg) !in bits
if(mod(size_of,8).eq.0) then
size_of=size_of/8 !in bytes
else
size_of=-1
endif
return
end function size_of
!----------------------------------
subroutine small_ascii(STR)
!Lowercases all capital Enlgish letters in string <STR>
implicit none
character(*), intent(inout):: STR
integer:: i,j,k,l,m,n,ia1,ia2,ia3
ia1=iachar('A')
ia2=iachar('Z')
ia3=iachar('a')
l=len_trim(STR)
do k=1,l