C
C Program SGFMAT to generate the system of linear equations for C the complex-valued coefficients of the structural Gabor functions C in decomposing a given gridded real-valued quantity C C Version: 7.40 C Date: 2017, May 19 C C Coded by: Ludek Klimes C Department of Geophysics, Charles University Prague, C Ke Karlovu 3, 121 16 Praha 2, Czech Republic, C http://sw3d.cz/staff/klimes.htm C C....................................................................... C C Two unknown real-valued coefficients stand in for the complex-valued C coefficient of each odd Gabor function. The odd unknown real-valued C coefficient is the real part, and the even unknown real-valued C coefficient is the imaginary part. These two real-valued coefficients C simultaneously represent the complex-valued coefficient of the C successive even Gabor function because the two Gabor functions and C the corresponding complex-valued coefficients are complex-conjugate. C Each odd unknown real-valued coefficient thus corresponds to twice the C real part of the Gabor function of the same odd index, and each even C unknown real-valued coefficient corresponds to twice the imaginary C part of the Gabor function of the same even index. C C Program SGFMAT generates the sparse symmetric real-valued matrix of C the scalar products between the above mentioned real-valued functions C corresponding to individual unknown real-valued coefficients, C and the real-valued vector of the scalar products of these real-valued C functions with the given gridded real-valued quantity. C C....................................................................... C C Description of data files: C C Input data read from the standard input device (*): C The data are read by the list directed input (free format) and C consist of a single string 'SEP': C 'SEP'...String in apostrophes containing the name of the input C SEP parameter or history file with the input data. C No default, 'SEP' must be specified and cannot be blank. C C C Input data file 'SEP': C File 'SEP' has the form of the SEP C parameter file. The parameters, which do not differ from their C defaults, need not be specified in file 'SEP'. C Data specifying input and output files: C SGF='string'... Name of the input file with structural Gabor C functions. Input parameters N1, N2, N3 determine whether C the structural Gabor functions are 1-D, 2-D or 3-D. C The wavemumber components and matrix elements C corresponding to the direction in which the number of C gridpoints is Ni=1 should equal zero. C Description of file SGF. C Default: SGF='sgf.out' C GRD='string'... Name of the formatted input file with a gridded C real-valued quantity. The file contains N1*N2*N3 values. C The file is not read if SGFRHS=' '. C For general description of files with gridded data refer C to file forms.htm. C No default: If SGFRHS is not blank, GRD must be specified C and cannot be blank. C SGFMAT='string'... Name of the header file of the output sparse C symmetric square real-valued matrix of the scalar products C between the real-valued functions corresponding to C individual unknown real-valued coefficients. This matrix C represents the coefficients of the linear equations for C the unknown complex-valued coefficients. C If SGFMAT is blank, the file is not created. C For general description of the files with matrices refer C to file forms.htm. C Default: SGFMAT='sgfmat.out' C SGFRHS='string'... Name of the header file of the output C real-valued vector of the scalar products of the given C gridded real-valued quantity with the real-valued C functions corresponding to individual unknown real-valued C coefficients. This vector represents the right-hand sides C of the linear equations for the unknown complex-valued C coefficients. C If SGFRHS is blank, the file is not created. C For general description of the files with matrices refer C to file forms.htm. C Default: SGFRHS='sgfrhs.out' C Data specifying dimensions of the input grid: C N1,N2,N3=integers... Numbers of gridpoints along the X1,X2,X3 C axes, respectively. These numbers also determine whether C the structural Gabor functions are 1-D, 2-D or 3-D. C Defaults: N1=1, N2=1, N3=1 C O1,O2,O3=reals... Coordinates of the origin of the grid, i.e., C of the first gridpoint. C Defaults: O1=0, O2=0, O3=0 C D1,D2,D3=reals... Grid intervals along the X1,X2,X3 axes, C respectively. C Defaults: D1=1, D2=1, D3=1 C Numerical parameters: C NULL=positive real... Scalar product smaller than NULL times the C product of the norms of multiplied functions is deemed to C be zero. C Default: 0.000001 C RELAMP=positive real... Relative decay of the Gaussian envelope C at which the loop over the points of the input grid is C terminated. C The relative error due to this economizing roughly C corresponds to the value of RELAMP. C Default: RELAMP=0.001 C Form of the output files with matrices: C FORMM='string' ... Form of the output files with matrices. Allowed C values are FORMM='formatted' and FORMM='unformatted'. C Default: FORMM='formatted' C Optional parameters specifying the form of the quantities C written to the output formatted files: C MINDIG,MAXDIG=positive integers ... See the description in file C forms.for. C NUMLINM=positive integer ... Number of reals or integers in one C line of the output file. See the description in file C mat.for. C Value of undefined quantities: C UNDEF=real... The value to be used for undefined real quantities. C Default: UNDEF=undefined value used in forms.for C C======================================================================= C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working array: INTEGER IRAM(MRAM),NRAM EQUIVALENCE (IRAM,RAM) C C....................................................................... C C External functions and subroutines: EXTERNAL ERROR,RSEP1,RSEP3T,RSEP3I,RSEP3R,RARRAY,WMATH,WMATD EXTERNAL STORE C C Constants: REAL COEF1D,COEF2D,COEF3D PARAMETER (COEF1D=1.772453851) PARAMETER (COEF2D=3.141592654) PARAMETER (COEF3D=5.568327997) C C Filenames and parameters: CHARACTER*80 FSEP,FSGF,FGRD,FMAT,FRHS,FDAT INTEGER LU1 PARAMETER (LU1=1) C C Input data: INTEGER N1,N2,N3,N1N2N3,NQ,NZ,NSGF,NQNSGF,NDIM REAL O1,O2,O3,D1,D2,D3,ANULL,ANULOG,RELAMP,RELLOG,ZERO(12) CHARACTER*1 TEXT C C Output data: INTEGER IMAT,NELEM0,NELEM CHARACTER*3 SPARSE CHARACTER*4 SYMM CHARACTER*11 FORM C C Gabor function b (beta): REAL BX1,BX2,BX3,BK1,BK2,BK3 REAL BY11,BY12,BY22,BY13,BY23,BY33,BYDET REAL BR11,BR12,BR22,BR13,BR23,BR33 REAL BXX1,BXX2,BXX3,BKK1,BKK2,BKK3 C Gabor function a (alpha), or 2-D projection of 3-D Gabor function: REAL AX1,AX2,AK1,AK2,AK3 REAL AY11,AY12,AY22,AY13,AY23,AY33,AYDET REAL AR11,AR12,AR22,AR13,AR23,AR33 C Differences and sums (Delta x, Delta k, matrices Y and R): REAL DX1,DX2,DX3,DK1,DK2,DK3 REAL YR11,YR21,YR31,YR12,YR22,YR32,YR13,YR23,YR33 REAL YI11,YI12,YI22,YI13,YI23,YI33 REAL RR11,RR12,RR22,RR13,RR23,RR33 REAL RI11,RI12,RI22,RI13,RI23,RI33 REAL ABKK1,ABKK2,ABKK3 C Determinant of matrix Y: REAL YDETR,YDETI,YDET2 C Inverse matrix to Y multiplied by the determinant of Y: REAL ZR11,ZR12,ZR22,ZR13,ZR23,ZR33 REAL ZI11,ZI12,ZI22,ZI13,ZI23,ZI33 C Products of matrices and vectors: REAL RDXR1,RDXR2,RDXR3,RDXI1,RDXI2,RDXI3 REAL ZRDXR1,ZRDXR2,ZRDXR3,ZRDXI1,ZRDXI2,ZRDXI3 REAL DXZDXR,DXZDXI C Square root of the determinant of matrix Y: REAL YDSQR,YDSQI,YDSQ2 C Calculation of scalar product (ga,gb): REAL EXPR,EXPI,C,S,GAGBR(2),GAGBI(2),VALUE,VALMIN C Calculation of scalar product with the given grid: REAL EXP0R,EXP0I,EXP1R,EXP1I,EXP1MR,EXP1MI,EXP2R,EXP2I REAL SUMR,SUMI INTEGER J1,J2,J3,K1,K2,K3,M1 C C Auxiliary variables: INTEGER ISGF1,ISGF2,IQ1,IQ2 INTEGER I1,I2,I3,I REAL DET,AUX C C....................................................................... C C Reading main input data: WRITE(*,'(A)') '+SGFMAT: Enter input filename: ' FSEP=' ' READ (*,*) FSEP IF(FSEP.EQ.' ') THEN C SGFMAT-01 CALL ERROR('SGFMAT-01: No input file specified') C Input file in the form of the SEP (Stanford Exploration Project) C parameter or history file must be specified. C There is no default filename. END IF CALL RSEP1(LU1,FSEP) WRITE(*,'(A)') '+SGFMAT: Working... ' C C Reading input and output filenames: CALL RSEP3T('SGF' ,FSGF,'sgf.out') CALL RSEP3T('GRD' ,FGRD,' ') CALL RSEP3T('SGFMAT',FMAT,'sgfmat.out') CALL RSEP3T('SGFRHS',FRHS,'sgfrhs.out') IF(FSGF.EQ.' ') THEN C SGFMAT-02 CALL ERROR('SGFMAT-02: Blank name of input file SGF') END IF IF(FGRD.EQ.' '.AND.FRHS.NE.' ') THEN C SGFMAT-03 CALL ERROR('SGFMAT-03: Name of input file GRD not specified') C If fileneme SGFRHS is not blank, the name of file GRD must be C specified and cannot be blank. END IF C C Reading other input SEP parameters: CALL RSEP3I('N1',N1,1) CALL RSEP3I('N2',N2,1) CALL RSEP3I('N3',N3,1) CALL RSEP3R('O1',O1,0.0) CALL RSEP3R('O2',O2,0.0) CALL RSEP3R('O3',O3,0.0) CALL RSEP3R('D1',D1,1.0) CALL RSEP3R('D2',D2,1.0) CALL RSEP3R('D3',D3,1.0) N1N2N3=N1*N2*N3 CALL RSEP3R('NULL',ANULL,0.000001) ANULOG=-ALOG(ANULL) CALL RSEP3R('RELAMP',RELAMP,0.001) RELLOG=-ALOG(RELAMP) C C Determination of NDIM: C 1-D: NDIM=1,2,3 C 2-D: NDIM=4,5,6 C 3-D: NDIM=7 NDIM=-1 IF(N1.GT.1) NDIM=NDIM+2 IF(N2.GT.1) NDIM=NDIM+3 IF(N3.GT.1) NDIM=NDIM+4 IF(NDIM.EQ.-1) THEN C SGFMAT-04 CALL ERROR('SGFMAT-04: N1=N2=N3=1 is not allowed') END IF IF(NDIM.EQ.2) THEN N1=N2 N2=1 O1=O2 D1=D2 ELSE IF(NDIM.EQ.3) THEN N1=N3 N3=1 O1=O3 D1=D3 ELSE IF(NDIM.EQ.5) THEN N2=N3 N3=1 O2=O3 D2=D3 ELSE IF(NDIM.EQ.6) THEN N1=N2 N2=N3 N3=1 O1=O2 O2=O3 D1=D2 D2=D3 ELSE IF(NDIM.EQ.8) THEN NDIM=7 END IF C C Determination of NQ and NZ: C NQ is the number of reals stored for each Gabor function C NZ is the number of input zeros for each Gabor function IF(NDIM.LE.3) THEN C 1-D NQ=6 NZ=12 ELSE IF(NDIM.LE.6) THEN C 2-D NQ=14 NZ=7 ELSE C 3-D: NQ=24 NZ=0 END IF C C Reading the Gabor functions: OPEN(LU1,FILE=FSGF,FORM='FORMATTED',STATUS='OLD') READ(LU1,*) (TEXT,I=1,20) NQNSGF=0 10 CONTINUE IF(NQNSGF+NQ.GT.MRAM) THEN C SGFMAT-05 CALL ERROR * ('SGFMAT-05: Too small array RAM for Gabor functions') C The input parameters of Gabor functions do not fit into array C RAM(MRAM). C The number of parameters of each Gabor function is C 1-D: NQ=6 C 2-D: NQ=14 C 3-D: NQ=24 END IF TEXT='$' GO TO (11,12,13,14,15,16,17) NDIM C SGFMAT-06 CALL ERROR('SGFMAT-06: Wrong value of NDIM') 11 CONTINUE READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),AUX,AUX, * RAM(NQNSGF+2),ZERO(1),ZERO(2), * RAM(NQNSGF+3),(ZERO(I),I=3,7), * RAM(NQNSGF+4),(ZERO(I),I=8,12) GO TO 18 12 CONTINUE READ(LU1,*,END=20) TEXT,AUX,RAM(NQNSGF+1),AUX, * ZERO(1),RAM(NQNSGF+2),(ZERO(I),I=2,4), * RAM(NQNSGF+3),(ZERO(I),I=5,9), * RAM(NQNSGF+4),(ZERO(I),I=10,12) GO TO 18 13 CONTINUE READ(LU1,*,END=20) TEXT,AUX,AUX,RAM(NQNSGF+1), * ZERO(1),ZERO(2),RAM(NQNSGF+2),(ZERO(I),I=3,7), * RAM(NQNSGF+3),(ZERO(I),I=8,12), * RAM(NQNSGF+4) GO TO 18 14 CONTINUE READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),RAM(NQNSGF+2),AUX, * RAM(NQNSGF+3),RAM(NQNSGF+4),ZERO(1), * RAM(NQNSGF+5),RAM(NQNSGF+6),RAM(NQNSGF+7), * (ZERO(I),I=2,4), * RAM(NQNSGF+8),RAM(NQNSGF+9),RAM(NQNSGF+10), * (ZERO(I),I=5,7) GO TO 18 15 CONTINUE READ(LU1,*,END=20) TEXT,RAM(NQNSGF+1),AUX,RAM(NQNSGF+2), * RAM(NQNSGF+3),ZERO(1),RAM(NQNSGF+4), * RAM(NQNSGF+5),ZERO(2),ZERO(3), * RAM(NQNSGF+6),ZERO(4),RAM(NQNSGF+7), * RAM(NQNSGF+8),ZERO(5),ZERO(6), * RAM(NQNSGF+9),ZERO(7),RAM(NQNSGF+10) GO TO 18 16 CONTINUE READ(LU1,*,END=20) TEXT,AUX,RAM(NQNSGF+1),RAM(NQNSGF+2), * ZERO(1),RAM(NQNSGF+3),RAM(NQNSGF+4), * ZERO(2),ZERO(3),RAM(NQNSGF+5), * ZERO(4),RAM(NQNSGF+6),RAM(NQNSGF+7), * ZERO(5),ZERO(6),RAM(NQNSGF+8), * ZERO(7),RAM(NQNSGF+9),RAM(NQNSGF+10) GO TO 18 17 CONTINUE READ(LU1,*,END=20) TEXT,(RAM(I),I=NQNSGF+1,NQNSGF+18) GO TO 18 18 CONTINUE IF(TEXT.EQ.'$') GO TO 20 NQNSGF=NQNSGF+NQ DO 19 I=1,NZ IF(ZERO(I).NE.0.0) THEN C SGFMAT-07 CALL ERROR ('SGFMAT-07: Non-zero input quantity') C The input wavemumber components and matrix elements C corresponding to the direction in which the number of C gridpoints is Ni=1 should equal zero. END IF 19 CONTINUE GO TO 10 20 CONTINUE CLOSE(LU1) NSGF=NQNSGF/NQ C NSGF is the number of odd Gabor functions. C RAM(1:NQNSGF) contain the parameters of odd Gabor functions. C C Checking the positive definiteness of the real part of matrix K, C calculating the squares of maximum differences between coordinates C and wavenumber components, and halving the components of matrix K: IF(NDIM.LE.3) THEN C 1-D: DO 32 I=0,NQNSGF-NQ,NQ IF(RAM(I+3).LE.0.0) THEN C SGFMAT-08 CALL ERROR('SGFMAT-08: Indefinite real part of matrix K') END IF C Squares of maximum phase-space coordinate differences RAM(I+5)=2.*ANULOG/RAM(I+3) RAM(I+6)=2.*ANULOG*(RAM(I+3)+RAM(I+4)**2/RAM(I+3)) C Halving the components of matrix K DO 31 I1=I+3,I+4 RAM(I1)=0.5*RAM(I1) 31 CONTINUE 32 CONTINUE ELSE IF(NDIM.LE.6) THEN C 2-D: DO 42 I=0,NQNSGF-NQ,NQ DET=RAM(I+5)*RAM(I+7)-RAM(I+6)**2 IF(RAM(I+5).LE.0.0.OR.DET.LE.0.0) THEN C SGFMAT-09 CALL ERROR('SGFMAT-09: Indefinite real part of matrix K') END IF C Squares of maximum phase-space coordinate differences RAM(I+11)=2.*ANULOG*RAM(I+7)/DET RAM(I+12)=2.*ANULOG*RAM(I+5)/DET AUX=RAM(I+8)*RAM(I+7)*RAM(I+8)+RAM(I+9)*RAM(I+5)*RAM(I+9) * -2.0*RAM(I+8)*RAM(I+6)*RAM(I+9) RAM(I+13)=2.*ANULOG*(RAM(I+5)+AUX/DET) AUX=RAM(I+9)*RAM(I+7)*RAM(I+9)+RAM(I+10)*RAM(I+5)*RAM(I+10) * -2.0*RAM(I+9)*RAM(I+6)*RAM(I+10) RAM(I+14)=2.*ANULOG*(RAM(I+7)+AUX/DET) C Halving the components of matrix K DO 41 I1=I+5,I+10 RAM(I1)=0.5*RAM(I1) 41 CONTINUE 42 CONTINUE ELSE C 3-D: DO 52 I=0,NQNSGF-NQ,NQ C Matrix Y BY11=RAM(I+7) BY12=RAM(I+8) BY22=RAM(I+9) BY13=RAM(I+10) BY23=RAM(I+11) BY33=RAM(I+12) C Matrix R BR11=-RAM(I+13) BR12=-RAM(I+14) BR22=-RAM(I+15) BR13=-RAM(I+16) BR23=-RAM(I+17) BR33=-RAM(I+18) C Inverse matrix to Y multiplied by det(Y) AY11=BY22*BY33-BY23*BY23 AY12=BY13*BY23-BY12*BY33 AY22=BY11*BY33-BY13*BY13 AY13=BY12*BY23-BY13*BY22 AY23=BY12*BY13-BY23*BY11 AY33=BY11*BY22-BY12*BY12 C det(Y) BYDET=BY11*AY11+BY12*AY12+BY13*AY13 IF(BY11.LE.0.0.OR.AY33.LE.0.0.OR.BYDET.LE.0.0) THEN C SGFMAT-10 CALL ERROR('SGFMAT-10: Indefinite real part of matrix K') END IF C Products of matrices YR11=AY11*BR11+AY12*BR12+AY13*BR13 YR21=AY12*BR11+AY22*BR12+AY23*BR13 YR31=AY13*BR11+AY23*BR12+AY33*BR13 YR12=AY11*BR12+AY12*BR22+AY13*BR23 YR22=AY12*BR12+AY22*BR22+AY23*BR23 YR32=AY13*BR12+AY23*BR22+AY33*BR23 YR13=AY11*BR12+AY12*BR22+AY13*BR23 YR23=AY12*BR13+AY22*BR23+AY23*BR33 YR33=AY13*BR13+AY23*BR23+AY33*BR33 AR11=BR11*YR11+BR12*YR21+BR13*YR31 AR22=BR12*YR12+BR22*YR22+BR23*YR32 AR33=BR13*YR13+BR23*YR23+BR33*YR33 C Squares of maximum phase-space coordinate differences RAM(I+19)=2.*ANULOG*AY11/BYDET RAM(I+20)=2.*ANULOG*AY22/BYDET RAM(I+21)=2.*ANULOG*AY33/BYDET RAM(I+22)=2.*ANULOG*(BY11+AR11/BYDET) RAM(I+23)=2.*ANULOG*(BY22+AR22/BYDET) RAM(I+24)=2.*ANULOG*(BY33+AR33/BYDET) C Halving the components of matrix K DO 51 I1=I+7,I+18 RAM(I1)=0.5*RAM(I1) 51 CONTINUE 52 CONTINUE END IF C C Check for the memory required for the calculation of the vector: IF(FRHS.NE.' '.AND.NQNSGF+N1N2N3.GT.MRAM) THEN C SGFMAT-11 CALL ERROR('SGFMAT-11: Too small array RAM for the grid') C The input parameters of Gabor functions and the input grid do C not fit together into array C RAM(MRAM). C The number of grid values is N1*N2*N3. C The number of parameters of each Gabor function is C 1-D: NQ=6 C 2-D: NQ=14 C 3-D: NQ=24 END IF C C....................................................................... C C Calculation of the symmetric matrix: IF(FMAT.NE.' ') THEN C C Preparation for storing the sparse matrix: NRAM=MRAM IMAT=NQNSGF+1 NELEM0=IMAT+2*NSGF NELEM=NELEM0 IRAM(IMAT)=NELEM0+1 DO 60 I=IMAT+1,NELEM IRAM(I)=0 60 CONTINUE C C Loops over Gabor functions: IF(NDIM.LE.3) THEN C 1-D: DO 112 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7),A,I3,A)') * '+SGFMAT: Matrix column',2*ISGF2,' of',2*NSGF, * ' (',NINT(100.*FLOAT(NELEM-NELEM0)/FLOAT(NRAM-NELEM0)), * '% RAM)' END IF IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BK1=RAM(IQ2+2) C Half matrix Yb BY11=RAM(IQ2+3) BYDET=BY11 C Half matrix Rb BR11=-RAM(IQ2+4) BXX1=RAM(IQ2+5) BKK1=RAM(IQ2+6) DO 111 ISGF1=0,ISGF2 IQ1=NQ*ISGF1 DX1=BX1-RAM(IQ1+1) IF(DX1*DX1.LE.RAM(IQ1+5)+BXX1) THEN AK1=RAM(IQ1+2) DK1=ABS(BK1)-ABS(AK1) ABKK1=RAM(IQ1+6)+BKK1 IF(DK1*DK1.LE.ABKK1) THEN GAGBR(1)=0.0 GAGBI(1)=0.0 GAGBR(2)=0.0 GAGBI(2)=0.0 C Half matrix Ya AY11=RAM(IQ1+3) AYDET=AY11 C Half matrix Ra AR11=-RAM(IQ1+4) C Real part of matrix Y YR11=BY11+AY11 C Imaginary part of matrix R RI11=BY11-AY11 DO 110 I=1,2 IF(I.EQ.2) THEN AK1=-AK1 END IF DK1=BK1-AK1 IF(DK1*DK1.LE.ABKK1) THEN IF(I.EQ.2) THEN AR11=-AR11 END IF C Imaginary part of matrix Y YI11=AR11-BR11 C Real part of matrix R RR11=AR11+BR11 C Determinant of matrix Y YDETR=YR11 YDETI=YI11 YDET2=YDETR*YDETR+YDETI*YDETI C Products of matrices and vectors RDXR1=RR11*DX1-DK1 RDXI1=RI11*DX1 DXZDXR=RDXR1*RDXR1-RDXI1*RDXI1 DXZDXI=RDXR1*RDXI1+RDXI1*RDXR1 C Real part of the exponent of scalar product EXPR=(DXZDXR*YDETR+DXZDXI*YDETI)/YDET2 EXPR=EXPR+DX1*YR11*DX1 EXPR=0.25*EXPR IF(EXPR.LE.ANULOG) THEN C Imaginary part of the exponent of scalar product EXPI=(DXZDXI*YDETR-DXZDXR*YDETI)/YDET2 EXPI=EXPI+DX1*YI11*DX1 EXPI=0.25*EXPI EXPI=EXPI+0.5*(BK1+AK1)*DX1 C=COS(-EXPI) S=SIN(-EXPI) C Square root of the determinant of matrix Y C (YDSQR^2-YDSQI^2=YDETR, 2*YDSQR*YDSQI=YDETI) YDSQ2=SQRT(YDET2) YDSQR=SQRT(0.5*(YDSQ2+YDETR)) YDSQI=SQRT(0.5*(YDSQ2-YDETR)) YDSQI=SIGN(YDSQI,YDETI) C Calculation of scalar product (ga,gb) AUX=EXP(-EXPR)/YDSQ2 GAGBR(I)=(YDSQR*C+YDSQI*S)*AUX GAGBI(I)=(YDSQR*S-YDSQI*C)*AUX END IF END IF 110 CONTINUE C C Storing the elements of the sparse matrix: C (ra-iya,rb+iyb)=(ra,rb)+(ya,yb)+i[(ra,yb)-(ya,rb)] C (ra+iya,rb+iyb)=(ra,rb)-(ya,yb)+i[(ra,yb)+(ya,rb)] AUX=2.0*SQRT(2.0*SQRT(AYDET*BYDET)) VALMIN=2.0*ANULL C Storing (2ra,2rb) I1=2*ISGF1+1 I2=2*ISGF2+1 VALUE=AUX*(GAGBR(1)+GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,2rb) I1=I1+1 IF(ISGF1.NE.ISGF2) THEN VALUE=AUX*(GAGBI(1)-GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF C Storing (2ra,-2yb) I1=I1-1 I2=I2+1 VALUE=-AUX*(GAGBI(1)+GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,-2yb) I1=I1+1 VALUE=AUX*(GAGBR(1)-GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF END IF 111 CONTINUE 112 CONTINUE ELSE IF(NDIM.LE.6) THEN C 2-D: DO 122 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7),A,I3,A)') * '+SGFMAT: Matrix column',2*ISGF2,' of',2*NSGF, * ' (',NINT(100.*FLOAT(NELEM-NELEM0)/FLOAT(NRAM-NELEM0)), * '% RAM)' END IF IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BX2=RAM(IQ2+2) BK1=RAM(IQ2+3) BK2=RAM(IQ2+4) C Half matrix Yb BY11=RAM(IQ2+5) BY12=RAM(IQ2+6) BY22=RAM(IQ2+7) BYDET=BY11*BY22-BY12*BY12 C Half matrix Rb BR11=-RAM(IQ2+8) BR12=-RAM(IQ2+9) BR22=-RAM(IQ2+10) BXX1=RAM(IQ2+11) BXX2=RAM(IQ2+12) BKK1=RAM(IQ2+13) BKK2=RAM(IQ2+14) DO 121 ISGF1=0,ISGF2 IQ1=NQ*ISGF1 DX1=BX1-RAM(IQ1+1) IF(DX1*DX1.LE.RAM(IQ1+11)+BXX1) THEN DX2=BX2-RAM(IQ1+2) IF(DX2*DX2.LE.RAM(IQ1+12)+BXX2) THEN AK1=RAM(IQ1+3) AK2=RAM(IQ1+4) DK1=ABS(BK1)-ABS(AK1) ABKK1=RAM(IQ1+13)+BKK1 IF(DK1*DK1.LE.ABKK1) THEN DK2=ABS(BK2)-ABS(AK2) ABKK2=RAM(IQ1+14)+BKK2 IF(DK2*DK2.LE.ABKK2) THEN GAGBR(1)=0.0 GAGBI(1)=0.0 GAGBR(2)=0.0 GAGBI(2)=0.0 C Half matrix Ya AY11=RAM(IQ1+5) AY12=RAM(IQ1+6) AY22=RAM(IQ1+7) AYDET=AY11*AY22-AY12*AY12 C Half matrix Ra AR11=-RAM(IQ1+8) AR12=-RAM(IQ1+9) AR22=-RAM(IQ1+10) C Real part of matrix Y YR11=BY11+AY11 YR12=BY12+AY12 YR22=BY22+AY22 C Imaginary part of matrix R RI11=BY11-AY11 RI12=BY12-AY12 RI22=BY22-AY22 DO 120 I=1,2 IF(I.EQ.2) THEN AK1=-AK1 AK2=-AK2 END IF DK1=BK1-AK1 IF(DK1*DK1.LE.ABKK1) THEN DK2=BK2-AK2 IF(DK2*DK2.LE.ABKK2) THEN IF(I.EQ.2) THEN AR11=-AR11 AR12=-AR12 AR22=-AR22 END IF C Imaginary part of matrix Y YI11=AR11-BR11 YI12=AR12-BR12 YI22=AR22-BR22 C Real part of matrix R RR11=AR11+BR11 RR12=AR12+BR12 RR22=AR22+BR22 C Inverse matrix to Y multiplied by det(Y) ZR11=YR22 ZR12=-YR12 ZR22=YR11 ZI11=YI22 ZI12=-YI12 ZI22=YI11 C Determinant of matrix Y YDETR=YR11*YR22-YI11*YI22-YR12*YR12+YI12*YI12 YDETI=YR11*YI22+YI11*YR22-YR12*YI12-YI12*YR12 YDET2=YDETR*YDETR+YDETI*YDETI C Products of matrices and vectors RDXR1=RR11*DX1+RR12*DX2-DK1 RDXR2=RR12*DX1+RR22*DX2-DK2 RDXI1=RI11*DX1+RI12*DX2 RDXI2=RI12*DX1+RI22*DX2 ZRDXR1=ZR11*RDXR1+ZR12*RDXR2 * -ZI11*RDXI1-ZI12*RDXI2 ZRDXR2=ZR12*RDXR1+ZR22*RDXR2 * -ZI12*RDXI1-ZI22*RDXI2 ZRDXI1=ZR11*RDXI1+ZR12*RDXI2 * +ZI11*RDXR1+ZI12*RDXR2 ZRDXI2=ZR12*RDXI1+ZR22*RDXI2 * +ZI12*RDXR1+ZI22*RDXR2 DXZDXR=RDXR1*ZRDXR1+RDXR2*ZRDXR2 * -RDXI1*ZRDXI1-RDXI2*ZRDXI2 DXZDXI=RDXR1*ZRDXI1+RDXR2*ZRDXI2 * +RDXI1*ZRDXR1+RDXI2*ZRDXR2 C Real part of the exponent of scalar product EXPR=(DXZDXR*YDETR+DXZDXI*YDETI)/YDET2 EXPR=EXPR+DX1*(YR11*DX1+2.*YR12*DX2) EXPR=EXPR+DX2* YR22*DX2 EXPR=0.25*EXPR IF(EXPR.LE.ANULOG) THEN C Imaginary part of the exponent of scalar product EXPI=(DXZDXI*YDETR-DXZDXR*YDETI)/YDET2 EXPI=EXPI+DX1*(YI11*DX1+2.* YI12*DX2) EXPI=EXPI+DX2* YI22*DX2 EXPI=0.25*EXPI EXPI=EXPI+0.5*((BK1+AK1)*DX1+(BK2+AK2)*DX2) C=COS(-EXPI) S=SIN(-EXPI) C Square root of the determinant of matrix Y C (YDSQR^2-YDSQI^2=YDETR, 2*YDSQR*YDSQI=YDETI) YDSQ2=SQRT(YDET2) YDSQR=SQRT(0.5*(YDSQ2+YDETR)) YDSQI=SQRT(0.5*(YDSQ2-YDETR)) YDSQI=SIGN(YDSQI,YDETI) C Calculation of scalar product (ga,gb) AUX=EXP(-EXPR)/YDSQ2 GAGBR(I)=(YDSQR*C+YDSQI*S)*AUX GAGBI(I)=(YDSQR*S-YDSQI*C)*AUX END IF END IF END IF 120 CONTINUE C C Storing the elements of the sparse matrix: C (ra-iya,rb+iyb)=(ra,rb)+(ya,yb)+i[(ra,yb)-(ya,rb)] C (ra+iya,rb+iyb)=(ra,rb)-(ya,yb)+i[(ra,yb)+(ya,rb)] AUX=2.0*SQRT(4.0*SQRT(AYDET*BYDET)) VALMIN=2.0*ANULL C Storing (2ra,2rb) I1=2*ISGF1+1 I2=2*ISGF2+1 VALUE=AUX*(GAGBR(1)+GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,2rb) I1=I1+1 IF(ISGF1.NE.ISGF2) THEN VALUE=AUX*(GAGBI(1)-GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF C Storing (2ra,-2yb) I1=I1-1 I2=I2+1 VALUE=-AUX*(GAGBI(1)+GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,-2yb) I1=I1+1 VALUE=AUX*(GAGBR(1)-GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF END IF END IF END IF 121 CONTINUE 122 CONTINUE ELSE C 3-D: DO 132 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7),A,I3,A)') * '+SGFMAT: Matrix column',2*ISGF2,' of',2*NSGF, * ' (',NINT(100.*FLOAT(NELEM-NELEM0)/FLOAT(NRAM-NELEM0)), * '% RAM)' END IF IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BX2=RAM(IQ2+2) BX3=RAM(IQ2+3) BK1=RAM(IQ2+4) BK2=RAM(IQ2+5) BK3=RAM(IQ2+6) C Half matrix Yb BY11=RAM(IQ2+7) BY12=RAM(IQ2+8) BY22=RAM(IQ2+9) BY13=RAM(IQ2+10) BY23=RAM(IQ2+11) BY33=RAM(IQ2+12) BYDET=BY11*(BY22*BY33-BY23*BY23) * +BY12*(BY13*BY23-BY33*BY12) * +BY13*(BY12*BY23-BY22*BY13) C Half matrix Rb BR11=-RAM(IQ2+13) BR12=-RAM(IQ2+14) BR22=-RAM(IQ2+15) BR13=-RAM(IQ2+16) BR23=-RAM(IQ2+17) BR33=-RAM(IQ2+18) BXX1=RAM(IQ2+19) BXX2=RAM(IQ2+20) BXX3=RAM(IQ2+21) BKK1=RAM(IQ2+22) BKK2=RAM(IQ2+23) BKK3=RAM(IQ2+24) DO 131 ISGF1=0,ISGF2 IQ1=NQ*ISGF1 DX1=BX1-RAM(IQ1+1) IF(DX1*DX1.LE.RAM(IQ1+19)+BXX1) THEN DX2=BX2-RAM(IQ1+2) IF(DX2*DX2.LE.RAM(IQ1+20)+BXX2) THEN DX3=BX3-RAM(IQ1+3) IF(DX3*DX3.LE.RAM(IQ1+21)+BXX3) THEN AK1=RAM(IQ1+4) AK2=RAM(IQ1+5) AK3=RAM(IQ1+6) DK1=ABS(BK1)-ABS(AK1) ABKK1=RAM(IQ1+22)+BKK1 IF(DK1*DK1.LE.ABKK1) THEN DK2=ABS(BK2)-ABS(AK2) ABKK2=RAM(IQ1+23)+BKK2 IF(DK2*DK2.LE.ABKK2) THEN DK3=ABS(BK3)-ABS(AK3) ABKK3=RAM(IQ1+24)+BKK3 IF(DK3*DK3.LE.ABKK3) THEN GAGBR(1)=0.0 GAGBI(1)=0.0 GAGBR(2)=0.0 GAGBI(2)=0.0 C Half matrix Ya AY11=RAM(IQ1+7) AY12=RAM(IQ1+8) AY22=RAM(IQ1+9) AY13=RAM(IQ1+10) AY23=RAM(IQ1+11) AY33=RAM(IQ1+12) AYDET=AY11*(AY22*AY33-AY23*AY23) * +AY12*(AY13*AY23-AY33*AY12) * +AY13*(AY12*AY23-AY22*AY13) C Half matrix Ra AR11=-RAM(IQ1+13) AR12=-RAM(IQ1+14) AR22=-RAM(IQ1+15) AR13=-RAM(IQ1+16) AR23=-RAM(IQ1+17) AR33=-RAM(IQ1+18) C Real part of matrix Y YR11=BY11+AY11 YR12=BY12+AY12 YR22=BY22+AY22 YR13=BY13+AY13 YR23=BY23+AY23 YR33=BY33+AY33 C Imaginary part of matrix R RI11=BY11-AY11 RI12=BY12-AY12 RI22=BY22-AY22 RI13=BY13-AY13 RI23=BY23-AY23 RI33=BY33-AY33 DO 130 I=1,2 IF(I.EQ.2) THEN AK1=-AK1 AK2=-AK2 AK3=-AK3 END IF DK1=BK1-AK1 IF(DK1*DK1.LE.ABKK1) THEN DK2=BK2-AK2 IF(DK2*DK2.LE.ABKK2) THEN DK3=BK3-AK3 IF(DK3*DK3.LE.ABKK3) THEN IF(I.EQ.2) THEN AR11=-AR11 AR12=-AR12 AR22=-AR22 AR13=-AR13 AR23=-AR23 AR33=-AR33 END IF C Imaginary part of matrix Y YI11=AR11-BR11 YI12=AR12-BR12 YI22=AR22-BR22 YI13=AR13-BR13 YI23=AR23-BR23 YI33=AR33-BR33 C Real part of matrix R RR11=AR11+BR11 RR12=AR12+BR12 RR22=AR22+BR22 RR13=AR13+BR13 RR23=AR23+BR23 RR33=AR33+BR33 C Inverse matrix to Y multiplied by det(Y) ZR11=YR22*YR33-YR23*YR23-YI22*YI33+YI23*YI23 ZR12=YR13*YR23-YR12*YR33-YI13*YI23+YI12*YI33 ZR22=YR11*YR33-YR13*YR13-YI11*YI33+YI13*YI13 ZR13=YR12*YR23-YR13*YR22-YI12*YI23+YI13*YI22 ZR23=YR12*YR13-YR23*YR11-YI12*YI13+YI23*YI11 ZR33=YR11*YR22-YR12*YR12-YI11*YI22+YI12*YI12 ZI11=YR22*YI33-YR23*YI23+YI22*YR33-YI23*YR23 ZI12=YR13*YI23-YR12*YI33+YI13*YR23-YI12*YR33 ZI22=YR11*YI33-YR13*YI13+YI11*YR33-YI13*YR13 ZI13=YR12*YI23-YR13*YI22+YI12*YR23-YI13*YR22 ZI23=YR12*YI13-YR23*YI11+YI12*YR13-YI23*YR11 ZI33=YR11*YI22-YR12*YI12+YI11*YR22-YI12*YR12 C Determinant of matrix Y YDETR=ZR11*YR11+ZR12*YR12+ZR13*YR13 * -ZI11*YI11-ZI12*YI12-ZI13*YI13 YDETI=ZR11*YI11+ZR12*YI12+ZR13*YI13 * +ZI11*YR11+ZI12*YR12+ZI13*YR13 YDET2=YDETR*YDETR+YDETI*YDETI C Products of matrices and vectors RDXR1=RR11*DX1+RR12*DX2+RR13*DX3-DK1 RDXR2=RR12*DX1+RR22*DX2+RR23*DX3-DK2 RDXR3=RR13*DX1+RR23*DX2+RR33*DX3-DK3 RDXI1=RI11*DX1+RI12*DX2+RI13*DX3 RDXI2=RI12*DX1+RI22*DX2+RI23*DX3 RDXI3=RI13*DX1+RI23*DX2+RI33*DX3 ZRDXR1=ZR11*RDXR1+ZR12*RDXR2+ZR13*RDXR3 * -ZI11*RDXI1-ZI12*RDXI2-ZI13*RDXI3 ZRDXR2=ZR12*RDXR1+ZR22*RDXR2+ZR23*RDXR3 * -ZI12*RDXI1-ZI22*RDXI2-ZI23*RDXI3 ZRDXR3=ZR13*RDXR1+ZR23*RDXR2+ZR33*RDXR3 * -ZI13*RDXI1-ZI23*RDXI2-ZI33*RDXI3 ZRDXI1=ZR11*RDXI1+ZR12*RDXI2+ZR13*RDXI3 * +ZI11*RDXR1+ZI12*RDXR2+ZI13*RDXR3 ZRDXI2=ZR12*RDXI1+ZR22*RDXI2+ZR23*RDXI3 * +ZI12*RDXR1+ZI22*RDXR2+ZI23*RDXR3 ZRDXI3=ZR13*RDXI1+ZR23*RDXI2+ZR33*RDXI3 * +ZI13*RDXR1+ZI23*RDXR2+ZI33*RDXR3 DXZDXR=RDXR1*ZRDXR1+RDXR2*ZRDXR2+RDXR3*ZRDXR3 * -RDXI1*ZRDXI1-RDXI2*ZRDXI2-RDXI3*ZRDXI3 DXZDXI=RDXR1*ZRDXI1+RDXR2*ZRDXI2+RDXR3*ZRDXI3 * +RDXI1*ZRDXR1+RDXI2*ZRDXR2+RDXI3*ZRDXR3 C Real part of the exponent of scalar product EXPR=(DXZDXR*YDETR+DXZDXI*YDETI)/YDET2 EXPR=EXPR+DX1*(YR11*DX1+2.*(YR12*DX2+YR13*DX3)) EXPR=EXPR+DX2*(YR22*DX2+2.* YR23*DX3) EXPR=EXPR+DX3* YR33*DX3 EXPR=0.25*EXPR IF(EXPR.LE.ANULOG) THEN C Imaginary part of the exponent of scalar product EXPI=(DXZDXI*YDETR-DXZDXR*YDETI)/YDET2 EXPI=EXPI+DX1*(YI11*DX1+2.*(YI12*DX2+YI13*DX3)) EXPI=EXPI+DX2*(YI22*DX2+2.* YI23*DX3) EXPI=EXPI+DX3* YI33*DX3 EXPI=0.25*EXPI EXPI=EXPI+0.5*((BK1+AK1)*DX1+(BK2+AK2)*DX2 * +(BK3+AK3)*DX3) C=COS(-EXPI) S=SIN(-EXPI) C Square root of the determinant of matrix Y C (YDSQR^2-YDSQI^2=YDETR, 2*YDSQR*YDSQI=YDETI) YDSQ2=SQRT(YDET2) YDSQR=SQRT(0.5*(YDSQ2+YDETR)) YDSQI=SQRT(0.5*(YDSQ2-YDETR)) YDSQI=SIGN(YDSQI,YDETI) IF(YDETR.LE.0.0) THEN IF(YDSQI*(ZI11+ZI22+ZI33).LT.0.0) THEN YDSQR=-YDSQR YDSQI=-YDSQI END IF END IF C Calculation of scalar product (ga,gb) AUX=EXP(-EXPR)/YDSQ2 GAGBR(I)=(YDSQR*C+YDSQI*S)*AUX GAGBI(I)=(YDSQR*S-YDSQI*C)*AUX END IF END IF END IF END IF 130 CONTINUE C C Storing the elements of the sparse matrix: C (ra-iya,rb+iyb)=(ra,rb)+(ya,yb)+i[(ra,yb)-(ya,rb)] C (ra+iya,rb+iyb)=(ra,rb)-(ya,yb)+i[(ra,yb)+(ya,rb)] AUX=2.0*SQRT(8.0*SQRT(AYDET*BYDET)) VALMIN=2.0*ANULL C Storing (2ra,2rb) I1=2*ISGF1+1 I2=2*ISGF2+1 VALUE=AUX*(GAGBR(1)+GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,2rb) I1=I1+1 IF(ISGF1.NE.ISGF2) THEN VALUE=AUX*(GAGBI(1)-GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF C Storing (2ra,-2yb) I1=I1-1 I2=I2+1 VALUE=-AUX*(GAGBI(1)+GAGBI(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) C Storing (-2ya,-2yb) I1=I1+1 VALUE=AUX*(GAGBR(1)-GAGBR(2)) CALL STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) END IF END IF END IF END IF END IF END IF 131 CONTINUE 132 CONTINUE END IF NELEM=NELEM-NELEM0 NELEM=NELEM/2 C NELEM is the number of non-zero elements of the sparse matrix. C C Writing the output sparse matrix: WRITE(*,'(A)') '+SGFMAT: Writing... ' FDAT=' ' SPARSE='CSC' SYMM='sym' FORM=' ' CALL WMATH(LU1,FMAT,FDAT,2*NSGF,2*NSGF,SPARSE,NELEM,SYMM,FORM) CALL WMATD(LU1,FDAT,2*NSGF,2*NSGF,SPARSE,NELEM,FORM,IMAT) C END IF C C....................................................................... C C Calculation of the vector: IF(FRHS.NE.' ') THEN FORM='FORMATTED' CALL RARRAY(LU1,FGRD,FORM,.TRUE.,0.0,N1N2N3,RAM(NQNSGF+1)) C IF(NDIM.LE.3) THEN C 1-D: DO 214 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7))') * '+SGFMAT: Vector component',2*ISGF2,' of',2*NSGF END IF C Initalizing the integrals SUMR=0.0 SUMI=0.0 C Quantities describing the Gabor function IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BK1=RAM(IQ2+2) C Half matrix Y BY11=RAM(IQ2+3) C Half matrix R BR11=-RAM(IQ2+4) C Determinant of half matrix Y BYDET=BY11 C Extent of the Gabor function along axis X1 AUX=SQRT(RELLOG/BY11) J1=MAX0(INT((BX1-AUX)/D1+0.999),0) K1=MIN0(INT((BX1+AUX)/D1+0.001),N1-1) IF(J1.LE.K1) THEN C Index of the central point along the gridline M1=MAX0(J1,MIN0(NINT(BX1/D1),K1)) C Relative coordinate of the central point DX1=O1+D1*FLOAT(M1)-BX1 C Exponent at the central point EXP0R=DX1*BY11*DX1 EXP0I=DX1*(BK1+BR11*DX1) C The first derivative of the exponent EXP1R=BY11*DX1 EXP1I=BR11*DX1 EXP1R=D1*(EXP1R+EXP1R) EXP1I=D1*(EXP1I+EXP1I+BK1) C Half the second derivative of the exponent EXP2R=D1*BY11*D1 EXP2I=D1*BR11*D1 C Half the second derivative minus the first derivative EXP1MR=EXP2R-EXP1R EXP1MI=EXP2I-EXP1I C Half the second derivative plus the first derivative EXP1R=EXP2R+EXP1R EXP1I=EXP2I+EXP1I C Second derivative of the exponent EXP2R=EXP2R+EXP2R EXP2I=EXP2I+EXP2I C Exponential function at the central point AUX=EXP(-EXP0R) C=COS(EXP0I) S=SIN(EXP0I) EXP0R=C*AUX EXP0I=S*AUX C Exponential correction at the central point (+) AUX=EXP(-EXP1R) C=COS(EXP1I) S=SIN(EXP1I) EXP1R=C*AUX EXP1I=S*AUX C Exponential correction at the central point (-) AUX=EXP(-EXP1MR) C=COS(EXP1MI) S=SIN(EXP1MI) EXP1MR=C*AUX EXP1MI=S*AUX C Constant correction to the correction AUX=EXP(-EXP2R) C=COS(EXP2I) S=SIN(EXP2I) EXP2R=C*AUX EXP2I=S*AUX C C Contribution to the integral at the central point I=NQNSGF+1+M1 AUX=RAM(I) SUMR=SUMR+AUX*EXP0R SUMI=SUMI+AUX*EXP0I C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C* Test for debugging * WRITE(*,*) C* C Loop over the gridpoints DO 210 I1=I+1,I+K1-M1 C* Test for debugging * DX1=O1+D1*FLOAT(I1-I-1+M1)-BX1 * CHKR=DX1*BY11*DX1 * CHKI=DX1*(BK1+BR11*DX1) * AUX=EXP(-CHKR) * C=COS(CHKI) * S=SIN(CHKI) * CHKR=C*AUX * CHKI=S*AUX * WRITE(*,*) CHKR,EXPR,CHKI,EXPI C* AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 210 CONTINUE C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C Exponential correction at the central point EXP1R=EXP1MR EXP1I=EXP1MI C Loop over the gridpoints DO 211 I1=I-1,I+J1-M1,-1 AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 211 CONTINUE END IF AUX=2.0*D1*SQRT(SQRT(2.0*BYDET)/COEF1D) RAM(2*ISGF2+1)= AUX*SUMR RAM(2*ISGF2+2)=-AUX*SUMI 214 CONTINUE ELSE IF(NDIM.LE.6) THEN C 2-D: DO 224 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7))') * '+SGFMAT: Vector component',2*ISGF2,' of',2*NSGF END IF C Initalizing the integrals SUMR=0.0 SUMI=0.0 C Quantities describing the Gabor function IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BX2=RAM(IQ2+2) BK1=RAM(IQ2+3) BK2=RAM(IQ2+4) C Half matrix Y BY11=RAM(IQ2+5) BY12=RAM(IQ2+6) BY22=RAM(IQ2+7) C Half matrix R BR11=-RAM(IQ2+8) BR12=-RAM(IQ2+9) BR22=-RAM(IQ2+10) C Determinant of half matrix Y BYDET=BY11*BY22-BY12*BY12 C Extent of the Gabor function along axis X2 AUX=SQRT(RELLOG*BY11/BYDET) J2=MAX0(INT((BX2-AUX-O2)/D2+0.999),0) K2=MIN0(INT((BX2+AUX-O2)/D2+0.001),N2-1) DO 222 I2=J2,K2 DX2=O2+D2*FLOAT(I2)-BX2 C Extent of the Gabor function along axis X1 AUX=RELLOG-BY22*DX2*DX2 IF(AUX.GE.0.0) THEN C Halfwidth of the Gabor function along axis X1 AUX=SQRT(AUX/BY11) C Central point along the gridline DX1=BX1-BY12*DX2/BY11-O1 J1=MAX0(INT((DX1-AUX)/D1+0.999),0) K1=MIN0(INT((DX1+AUX)/D1+0.001),N1-1) IF(J1.LE.K1) THEN C Index of the central point along the gridline M1=MAX0(J1,MIN0(NINT(DX1/D1),K1)) C Relative coordinate of the central point DX1=O1+D1*FLOAT(M1)-BX1 C Exponent at the central point EXP0R=DX1*(BY11*DX1+2.0*BY12*DX2)+DX2*BY22*DX2 EXP0I= DX1*(BK1+BR11*DX1+2.0*BR12*DX2) EXP0I=EXP0I+DX2*(BK2+BR22*DX2) C The first derivative of the exponent EXP1R=BY11*DX1+BY12*DX2 EXP1I=BR11*DX1+BR12*DX2 EXP1R=D1*(EXP1R+EXP1R) EXP1I=D1*(EXP1I+EXP1I+BK1) C Half the second derivative of the exponent EXP2R=D1*BY11*D1 EXP2I=D1*BR11*D1 C Half the second derivative minus the first derivative EXP1MR=EXP2R-EXP1R EXP1MI=EXP2I-EXP1I C Half the second derivative plus the first derivative EXP1R=EXP2R+EXP1R EXP1I=EXP2I+EXP1I C Second derivative of the exponent EXP2R=EXP2R+EXP2R EXP2I=EXP2I+EXP2I C Exponential function at the central point AUX=EXP(-EXP0R) C=COS(EXP0I) S=SIN(EXP0I) EXP0R=C*AUX EXP0I=S*AUX C Exponential correction at the central point (+) AUX=EXP(-EXP1R) C=COS(EXP1I) S=SIN(EXP1I) EXP1R=C*AUX EXP1I=S*AUX C Exponential correction at the central point (-) AUX=EXP(-EXP1MR) C=COS(EXP1MI) S=SIN(EXP1MI) EXP1MR=C*AUX EXP1MI=S*AUX C Constant correction to the correction AUX=EXP(-EXP2R) C=COS(EXP2I) S=SIN(EXP2I) EXP2R=C*AUX EXP2I=S*AUX C C Contribution to the integral at the central point I=NQNSGF+1+M1+N1*I2 AUX=RAM(I) SUMR=SUMR+AUX*EXP0R SUMI=SUMI+AUX*EXP0I C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C* Test for debugging * WRITE(*,*) C* C Loop over the gridpoints DO 220 I1=I+1,I+K1-M1 C* Test for debugging * DX1=O1+D1*FLOAT(I1-I-1+M1)-BX1 * DX1=O1+D1*FLOAT(I1-I+M1)-BX1 * CHKR= DX1*(BY11*DX1+2.*BY12*DX2) * CHKR=CHKR+DX2* BY22*DX2 * CHKI= DX1*(BK1+BR11*DX1+2.*BR12*DX2) * CHKI=CHKI+DX2*(BK2+BR22*DX2) * AUX=EXP(-CHKR) * C=COS(CHKI) * S=SIN(CHKI) * CHKR=C*AUX * CHKI=S*AUX * WRITE(*,*) CHKR,EXPR,CHKI,EXPI * EXPR=CHKR * EXPI=CHKI C* AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 220 CONTINUE C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C Exponential correction at the central point EXP1R=EXP1MR EXP1I=EXP1MI C* Test for debugging * WRITE(*,*) C* C Loop over the gridpoints DO 221 I1=I-1,I+J1-M1,-1 C* Test for debugging * DX1=O1+D1*FLOAT(I1-I+1+M1)-BX1 * DX1=O1+D1*FLOAT(I1-I+M1)-BX1 * CHKR= DX1*(BY11*DX1+2.*BY12*DX2) * CHKR=CHKR+DX2* BY22*DX2 * CHKI= DX1*(BK1+BR11*DX1+2.*BR12*DX2) * CHKI=CHKI+DX2*(BK2+BR22*DX2) * AUX=EXP(-CHKR) * C=COS(CHKI) * S=SIN(CHKI) * CHKR=C*AUX * CHKI=S*AUX * WRITE(*,*) CHKR,EXPR,CHKI,EXPI * EXPR=CHKR * EXPI=CHKI C* AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 221 CONTINUE END IF END IF 222 CONTINUE AUX=2.0*D1*D2*SQRT(SQRT(4.0*BYDET)/COEF2D) RAM(2*ISGF2+1)= AUX*SUMR RAM(2*ISGF2+2)=-AUX*SUMI 224 CONTINUE ELSE C 3-D: DO 234 ISGF2=0,NSGF-1 IF(MOD(ISGF2,50).EQ.0) THEN WRITE(*,'(2(A,I7))') * '+SGFMAT: Vector component',2*ISGF2,' of',2*NSGF END IF C Initalizing the integrals SUMR=0.0 SUMI=0.0 C Quantities describing the Gabor function IQ2=NQ*ISGF2 BX1=RAM(IQ2+1) BX2=RAM(IQ2+2) BX3=RAM(IQ2+3) BK1=RAM(IQ2+4) BK2=RAM(IQ2+5) BK3=RAM(IQ2+6) C Half matrix Y BY11=RAM(IQ2+7) BY12=RAM(IQ2+8) BY22=RAM(IQ2+9) BY13=RAM(IQ2+10) BY23=RAM(IQ2+11) BY33=RAM(IQ2+12) C Half matrix R BR11=-RAM(IQ2+13) BR12=-RAM(IQ2+14) BR22=-RAM(IQ2+15) BR13=-RAM(IQ2+16) BR23=-RAM(IQ2+17) BR33=-RAM(IQ2+18) C Determinant of the 2*2 submatrix of half Y DET=BY11*BY22-BY12*BY12 C Matrix inverse to the 2*2 submatrix of Y times (BY13,BY23) AY13=( BY22*BY13-BY12*BY23)/DET AY23=(-BY12*BY13+BY11*BY23)/DET C Determinant of half matrix Y BYDET=BY33*DET-BY11*BY23*BY23 * +BY13*(2.0*BY12*BY23-BY22*BY13) C Extent of the Gabor function along axis X3 AUX=SQRT(RELLOG*DET/BYDET) J3=MAX0(INT((BX3-AUX-O3)/D3+0.999),0) K3=MIN0(INT((BX3+AUX-O3)/D3+0.001),N3-1) DO 233 I3=J3,K3 DX3=O3+D3*FLOAT(I3)-BX3 C C Transforming 3-D Gabor packet to 2-D Gabor packet C Shift of the central point DX1=-AY13*DX3 DX2=-AY23*DX3 C Central point of the 2-D Gabor packet AX1=BX1+DX1 AX2=BX2+DX2 C Wavenumber vetor of the 2-D Gabor packet AK1=BK1+2.0*(BR11*DX1+BR12*DX2+BR13*DX3) AK2=BK2+2.0*(BR12*DX1+BR22*DX2+BR23*DX3) C Exponent of the 2-D Gabor packet at its central point AR33= DX1*(BK1+BR11*DX1+2.0*(BR12*DX2+BR13*DX3)) AR33=AR33+DX2*(BK2+BR22*DX2+2.0* BR23*DX3) AR33=AR33+DX3*(BK3+BR33*DX3) AY33=(BY33-BY13*AY13-BY23*AY23)*DX3*DX3 C Matrices BY11,BY12,BY22 and BR11,BR12,BR22 keep unchanged. C C Extent of the Gabor function along axis X2 AUX=RELLOG-AY33 IF(AUX.GE.0.0) THEN AUX=SQRT(AUX*BY11/DET) J2=MAX0(INT((AX2-AUX-O2)/D2+0.999),0) K2=MIN0(INT((AX2+AUX-O2)/D2+0.001),N2-1) DO 232 I2=J2,K2 DX2=O2+D2*FLOAT(I2)-AX2 C Extent of the Gabor function along axis X1 AUX=RELLOG-AY33-(DET/BY11)*DX2*DX2 IF(AUX.GE.0.0) THEN AUX=SQRT(AUX/BY11) DX1=AX1-BY12*DX2/BY11-O1 J1=MAX0(INT((DX1-AUX)/D1+0.999),0) K1=MIN0(INT((DX1+AUX)/D1+0.001),N1-1) IF(J1.LE.K1) THEN C Index of the central point along the gridline M1=MAX0(J1,MIN0(NINT(DX1/D1),K1)) C Relative coordinate of the central point DX1=O1+D1*FLOAT(M1)-AX1 C Exponent at the central point EXP0R=AY33 EXP0R=EXP0R+DX1*(BY11*DX1+2.0*BY12*DX2) EXP0R=EXP0R+DX2* BY22*DX2 EXP0I=AR33 EXP0I=EXP0I+DX1*(AK1+BR11*DX1+2.0*BR12*DX2) EXP0I=EXP0I+DX2*(AK2+BR22*DX2) C The first derivative of the exponent EXP1R=BY11*DX1+BY12*DX2 EXP1I=BR11*DX1+BR12*DX2 EXP1R=D1*(EXP1R+EXP1R) EXP1I=D1*(EXP1I+EXP1I+AK1) C Half the second derivative of the exponent EXP2R=D1*BY11*D1 EXP2I=D1*BR11*D1 C Half the second derivative minus first derivative EXP1MR=EXP2R-EXP1R EXP1MI=EXP2I-EXP1I C Half the second derivative plus first derivative EXP1R=EXP2R+EXP1R EXP1I=EXP2I+EXP1I C Second derivative of the exponent EXP2R=EXP2R+EXP2R EXP2I=EXP2I+EXP2I C Exponential function at the central point AUX=EXP(-EXP0R) C=COS(EXP0I) S=SIN(EXP0I) EXP0R=C*AUX EXP0I=S*AUX C Exponential correction at the central point (+) AUX=EXP(-EXP1R) C=COS(EXP1I) S=SIN(EXP1I) EXP1R=C*AUX EXP1I=S*AUX C Exponential correction at the central point (-) AUX=EXP(-EXP1MR) C=COS(EXP1MI) S=SIN(EXP1MI) EXP1MR=C*AUX EXP1MI=S*AUX C Constant correction to the correction AUX=EXP(-EXP2R) C=COS(EXP2I) S=SIN(EXP2I) EXP2R=C*AUX EXP2I=S*AUX C C Contribution to the integral at the central point I=NQNSGF+1+M1+N1*(I2+N2*I3) AUX=RAM(I) SUMR=SUMR+AUX*EXP0R SUMI=SUMI+AUX*EXP0I C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C* Test for debugging * WRITE(*,*) * DX2=O2+D2*FLOAT(I2)-BX2 C* C Loop over the gridpoints DO 230 I1=I+1,I+K1-M1 C* Test for debugging * DX1=O1+D1*FLOAT(I1-I-1+M1)-BX1 * CHKR= DX1*(BY11*DX1+2.*(BY12*DX2+BY13*DX3)) * CHKR=CHKR+DX2*(BY22*DX2+2.* BY23*DX3) * CHKR=CHKR+DX3* BY33*DX3 * CHKI=BK1*DX1+BK2*DX2+BK3*DX3 * CHKI=CHKI+DX1*(BR11*DX1+2.*(BR12*DX2+BR13*DX3)) * CHKI=CHKI+DX2*(BR22*DX2+2.* BR23*DX3) * CHKI=CHKI+DX3* BR33*DX3 * AUX=EXP(-CHKR) * C=COS(CHKI) * S=SIN(CHKI) * CHKR=C*AUX * CHKI=S*AUX * WRITE(*,*) CHKR,EXPR,CHKI,EXPI C* AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 230 CONTINUE C C Exponential function at the central point EXPR=EXP0R EXPI=EXP0I C Exponential correction at the central point EXP1R=EXP1MR EXP1I=EXP1MI C* Test for debugging * WRITE(*,*) * DX2=O2+D2*FLOAT(I2)-BX2 C* C Loop over the gridpoints DO 231 I1=I-1,I+J1-M1,-1 C* Test for debugging * DX1=O1+D1*FLOAT(I1-I+1+M1)-BX1 * CHKR= DX1*(BY11*DX1+2.*(BY12*DX2+BY13*DX3)) * CHKR=CHKR+DX2*(BY22*DX2+2.* BY23*DX3) * CHKR=CHKR+DX3* BY33*DX3 * CHKI=BK1*DX1+BK2*DX2+BK3*DX3 * CHKI=CHKI+DX1*(BR11*DX1+2.*(BR12*DX2+BR13*DX3)) * CHKI=CHKI+DX2*(BR22*DX2+2.* BR23*DX3) * CHKI=CHKI+DX3* BR33*DX3 * AUX=EXP(-CHKR) * C=COS(CHKI) * S=SIN(CHKI) * CHKR=C*AUX * CHKI=S*AUX * WRITE(*,*) CHKR,EXPR,CHKI,EXPI C* AUX =EXPR*EXP1R-EXPI*EXP1I EXPI=EXPR*EXP1I+EXPI*EXP1R EXPR=AUX AUX =EXP1R*EXP2R-EXP1I*EXP2I EXP1I=EXP1R*EXP2I+EXP1I*EXP2R EXP1R=AUX AUX=RAM(I1) SUMR=SUMR+AUX*EXPR SUMI=SUMI+AUX*EXPI 231 CONTINUE END IF END IF 232 CONTINUE END IF 233 CONTINUE AUX=2.0*D1*D2*D3*SQRT(SQRT(8.0*BYDET)/COEF3D) RAM(2*ISGF2+1)= AUX*SUMR RAM(2*ISGF2+2)=-AUX*SUMI 234 CONTINUE END IF C C Writing the output vector: WRITE(*,'(A)') '+SGFMAT: Writing... ' FDAT=' ' SPARSE=' ' SYMM=' ' FORM=' ' CALL WMATH(LU1,FRHS,FDAT,2*NSGF,1,SPARSE,2*NSGF,SYMM,FORM) CALL WMATD(LU1,FDAT,2*NSGF,1,SPARSE,2*NSGF,FORM,1) C END IF C WRITE(*,'(A)') '+SGFMAT: Done. ' STOP END C C======================================================================= C C C SUBROUTINE STORE(I1,I2,VALUE,VALMIN,IMAT,NELEM,NRAM) INTEGER I1,I2,IMAT,NELEM,NRAM REAL VALUE,VALMIN C C Subroutine storing a given matrix element into arrays IRAM and RAM. C C Input: C I1... Row index of the element of a sparse matrix. C I2... Column index of the element of a sparse matrix. C VALUE.. Value of the element of a sparse matrix. C VALMIN. Minimim value for storing. C If ABS(VALUE).LT.VALMIN, nothing is stored. C IMAT... Address of the first storage location of the matrix C in array RAM of common block /RAMC/. C NELEM.. Index of the last stored value in array RAM. C NRAM... Dimension of the usable part of array RAM. C C Output if ABS(VALUE).GE.VALMIN: C NELEM.. Input value increased by 2. C Array RAM of common block /RAMC/ is updated. C C No output if ABS(VALUE).LT.VALMIN. C C----------------------------------------------------------------------- C C Common block /RAMC/: INCLUDE 'ram.inc' C ram.inc C C Allocation of working array: INTEGER IRAM(MRAM) EQUIVALENCE (IRAM,RAM) C C....................................................................... C C External functions and subroutines: EXTERNAL ERROR C C Auxiliary storage locations: INTEGER I,J C C----------------------------------------------------------------------- C IF(ABS(VALUE).GE.VALMIN) THEN IF(NELEM+2.GT.NRAM) THEN C SGFMAT-12 CALL ERROR('SGFMAT-12: Too small array RAM for the matrix') C The input parameters of Gabor functions and the elements C of the sparse symmetric matrix do not fit into array C RAM(MRAM). C The number of parameters of each Gabor function is C 1-D: NQ=6 C 2-D: NQ=14 C 3-D: NQ=24 C Dimension MRAM of array RAM should be at least NQ*NSF+3*NELEM, C where NSF is the number of odd Gabor functions and C NELEM is the number of non-zero elements of the sparse C symmetric matrix. C You may increase the value of input SEP parameter NULL, C increase dimension C MRAM and recompile C the program, or decrease the number of input Gabor functions. END IF NELEM=NELEM+2 C C Updating the length of the column I=IRAM(IMAT+I2) IF(I.EQ.0) THEN I=IRAM(IMAT+I2-1) IF(I.EQ.0) THEN I=IRAM(IMAT+I2-2) IF(I.EQ.0) THEN C SGFMAT-13 CALL ERROR('SGFMAT-13: Zero column of the matrix') C This error should not appear. Contact the authors. END IF IRAM(IMAT+I2-1)=I END IF END IF IRAM(IMAT+I2)=I+2 C C Moving the already stored elements of the next column J=IMAT+I2+1 IF(J.LT.IRAM(IMAT)) THEN J=IRAM(J) IF(J.NE.0) THEN IRAM(IMAT+I2+1)=J+2 DO 10 J=J,I+2,-2 IRAM(J )=IRAM(J-2) RAM (J+1)=RAM (J-1) 10 CONTINUE END IF END IF C C Storing the matrix element IRAM(I )=I1 RAM (I+1)=VALUE END IF RETURN END C C======================================================================= C INCLUDE 'error.for' C error.for INCLUDE 'sep.for' C sep.for INCLUDE 'length.for' C length.for INCLUDE 'forms.for' C forms.for INCLUDE 'mat.for' C mat.for C C======================================================================= C