C GRQUA3 SOURCE CHAT 05/01/13 00:21:31 5004 SUBROUTINE GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, & MATCEL,MLRSIG,MLRCEL,COEF,INDER) C----------------------------------------------------------------------- C PROJET : CASTEM 2000 C C NOM : GRQUA3 C C DESCRIPTION : Calcule les coefficients à appliquer C sur la fonction dont onv eut calculer le gradient et C le Hessien. C C LANGAGE : FORTRAN 77 + ESOPE 2000 C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF C C----------------------------------------------------------------------- C C APPELES (E/S) : C C C APPELES (Calcul): INVSVD C C----------------------------------------------------------------------- C C ENTREES : CONTOL = REAL 8 t.q. 1/cond(MATA) > CONTOL C Si non rank(MATA) < NVOIMI - 1 C C NBMAX, NTOTP, NVOIMI : INTEGERs qui interviennent C dans les dimensions de C certains matrices C C MLEVTO = pointeur sur un list d'entiers C (NBMAX - 1) * NBMAX : C MLEVTO.LECT(1) = point ou je veux calculer C gradient et hessien C MLEVTO.LECT(i), i=2, NTOTP = voisins C C MATA = pointeur d'une matrice C (NBMAXx(NBMAX-1),NVOIMI-1) qui contient la C matrice (NTOTP-1,NVOIMI-1), a pseudoinverser C C MATU = pointeur d'une matrice C (NBMAXx(NBMAX-1),NVOIMI-1) on met la matrice C (NTOTP-1,NVOIMI-1) des vectors singuliers C droites de MATA C C MATV = pointeur de la matrice (NVOIMI-1,NVOIMI-1) C des vectors singuliers gauches de MATA C C MATCEL = pointeur sur une matrice de travail C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient C momentanément b C la matrice (NTOTP-1,NVOIMI-1), transposé C de la pseudo-inverse de MATA C C MLRSIG = pointeur sur une liste de reels (NVOIMI-1) C qui contient lesvaleurs singuliers de MATA C C MLRCEL = pointeur sur une liste de reels C (NBMAXx(NBMAX-1)) (vecteur de travail) C C C SORTIES : COEF = pointeur sur une matrice C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient C la matrice (NTOTP,NVOIMI-1) pour calculer C le coefficients C C INDER = 1 -> probleme mal posé C C INDER = 2 -> erreur anormale C C----------------------------------------------------------------- C C HISTORIQUE (Anomalies et modifications eventuelles) C C HISTORIQUE : 10.05.2000, Creation C C----------------------------------------------------------------- C C IMPLICIT INTEGER(I-N) -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMLENTI -INC SMLREEL C C----------------------------------------------------------------- C POINTEUR MLEVTO.MLENTI POINTEUR MLRSIG.MLREEL, MLRCEL.MLREEL C SEGMENT MATRIX REAL*8 MAT(N1,N2) ENDSEGMENT C POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX, MATCEL.MATRIX, & COEF.MATRIX C C**** Variables de COOPTIO C C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES C & ,IECHO, IIMPI, IOSPI C & ,IDIM CC & ,MCOORD C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU C & ,NORINC,NORVAL,NORIND,NORVAD C & ,NUCROU, IPSAUV C INTEGER NBMAX, NTOTP, NVOIMI, INDER, NOCENE, ICOORX, NOCEN1, IVOIS & , ILIG, ICOL, IVECT, IERSVD REAL*8 CONTOL, HELE, XCEN, YCEN, XCEN1, YCEN1, DX, DY, DX2, DY2 & , HELE2, SVMAX, SVMIN C C INDER = 0 C NOCENE = MLEVTO.LECT(1) HELE = 0.0D0 ICOORX = ((IDIM + 1) * (NOCENE - 1))+1 XCEN = MCOORD.XCOOR(ICOORX) YCEN = MCOORD.XCOOR(ICOORX+1) C C**** La submatrice de MATA qui nous interesse est (NTOTP-1)x(NVOIMI-1) C Amelioration: ca serait mieux travailler avec sa transposé C Mais invsvd.eso a ete teste avec m > n C DO IVOIS = 1, NTOTP-1, 1 NOCEN1 = MLEVTO.LECT(IVOIS+1) ICOORX = ((IDIM + 1) * (NOCEN1 - 1))+1 XCEN1 = MCOORD.XCOOR(ICOORX) YCEN1 = MCOORD.XCOOR(ICOORX+1) DX = XCEN1 - XCEN DY = YCEN1 - YCEN DX2 = DX * DX DY2 = DY * DY HELE = MAX(HELE,SQRT(DX2+DY2)) MATA.MAT(IVOIS,1) = DX MATA.MAT(IVOIS,2) = DY MATA.MAT(IVOIS,3) = 0.5D0 * DX2 MATA.MAT(IVOIS,4) = 0.5D0 * DY2 MATA.MAT(IVOIS,5) = DX * DY ENDDO C C******* On doit adimensioner MATA C HELE2 = HELE * HELE DO IVOIS = 1, NTOTP - 1, 1 MATA.MAT(IVOIS,1) = MATA.MAT(IVOIS,1) / HELE MATA.MAT(IVOIS,2) = MATA.MAT(IVOIS,2) / HELE MATA.MAT(IVOIS,3) = MATA.MAT(IVOIS,3) / HELE2 MATA.MAT(IVOIS,4) = MATA.MAT(IVOIS,4) / HELE2 MATA.MAT(IVOIS,5) = MATA.MAT(IVOIS,5) / HELE2 ENDDO C CALL INVSVD( NBMAX * (NBMAX -1),NTOTP-1,NVOIMI-1 ,MATA.MAT, & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD, & MLRCEL.PROG) C IF(IERSVD .NE. 0)THEN C OPEN(1,FILE='SVD.PBM') C WRITE(1,*) NBMAX * (NBMAX -1) C WRITE(1,*) NTOTP-1 C WRITE(1,*) NVOIMI-1 C WRITE(1,*) C DO IVOIS = 1, NTOTP - 1, 1 C DO ILIG = 1, NVOIMI - 1 C WRITE(1,*) MATA.MAT(IVOIS,ILIG) C ENDDO C ENDDO C CLOSE(1) INDER = 2 GOTO 9999 ENDIF C SVMAX = 0.0D0 DO IVECT = 1, NVOIMI-1, 1 SVMAX = MAX(SVMAX,MLRSIG.PROG(IVECT)) ENDDO C SVMIN = SVMAX DO IVECT = 1, NVOIMI-1, 1 SVMIN = MIN(SVMIN,MLRSIG.PROG(IVECT)) ENDDO C IF((SVMIN / SVMAX) .LT. CONTOL)THEN DO IVOIS = 1 , NTOTP DO IVECT = 1, 5 COEF.MAT(IVOIS,IVECT) = 0.0D0 ENDDO ENDDO INDER = 1 GOTO 9999 ENDIF C C******* On doit calculer MATCEL (transpose de la pseudo inverse de MATA) C DO ICOL = 1, NVOIMI - 1, 1 DO ILIG = 1, NTOTP - 1, 1 MATCEL.MAT(ILIG,ICOL) = 0.0D0 ENDDO ENDDO C DO IVECT = 1, NVOIMI - 1, 1 DO ICOL = 1, NVOIMI - 1, 1 DO ILIG = 1, NTOTP - 1, 1 MATCEL.MAT(ILIG,ICOL) = MATCEL.MAT(ILIG,ICOL) + & (MATU.MAT(ILIG,IVECT) * MATV.MAT(ICOL,IVECT) & / MLRSIG.PROG(IVECT)) ENDDO ENDDO ENDDO C DO IVECT = 1, 5 COEF.MAT(1,IVECT) = 0.0D0 ENDDO C DO IVECT = 1, 2 DO IVOIS = 1, NTOTP - 1, 1 COEF.MAT(1+IVOIS,IVECT) = & MATCEL.MAT(IVOIS,IVECT) / HELE COEF.MAT(1,IVECT) = & COEF.MAT(1,IVECT) - & COEF.MAT(1+IVOIS,IVECT) ENDDO ENDDO C DO IVECT = 3, 5 DO IVOIS = 1, NTOTP - 1, 1 COEF.MAT(1+IVOIS,IVECT) = & MATCEL.MAT(IVOIS,IVECT) / HELE2 COEF.MAT(1,IVECT) = & COEF.MAT(1,IVECT) - COEF.MAT(1+IVOIS,IVECT) ENDDO ENDDO 9999 RETURN C END