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