1 : C GRQUA3 SOURCE CHAT 05/01/13 00:21:31 5004 2 : SUBROUTINE GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, 3 : & MATCEL,MLRSIG,MLRCEL,COEF,INDER) 4 : C----------------------------------------------------------------------- 5 : C PROJET : CASTEM 2000 6 : C 7 : C NOM : GRQUA3 8 : C 9 : C DESCRIPTION : Calcule les coefficients à appliquer 10 : C sur la fonction dont onv eut calculer le gradient et 11 : C le Hessien. 12 : C 13 : C LANGAGE : FORTRAN 77 + ESOPE 2000 14 : C 15 : C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF 16 : C 17 : C----------------------------------------------------------------------- 18 : C 19 : C APPELES (E/S) : 20 : C 21 : C 22 : C APPELES (Calcul): INVSVD 23 : C 24 : C----------------------------------------------------------------------- 25 : C 26 : C ENTREES : CONTOL = REAL 8 t.q. 1/cond(MATA) > CONTOL 27 : C Si non rank(MATA) < NVOIMI - 1 28 : C 29 : C NBMAX, NTOTP, NVOIMI : INTEGERs qui interviennent 30 : C dans les dimensions de 31 : C certains matrices 32 : C 33 : C MLEVTO = pointeur sur un list d'entiers 34 : C (NBMAX - 1) * NBMAX : 35 : C MLEVTO.LECT(1) = point ou je veux calculer 36 : C gradient et hessien 37 : C MLEVTO.LECT(i), i=2, NTOTP = voisins 38 : C 39 : C MATA = pointeur d'une matrice 40 : C (NBMAXx(NBMAX-1),NVOIMI-1) qui contient la 41 : C matrice (NTOTP-1,NVOIMI-1), a pseudoinverser 42 : C 43 : C MATU = pointeur d'une matrice 44 : C (NBMAXx(NBMAX-1),NVOIMI-1) on met la matrice 45 : C (NTOTP-1,NVOIMI-1) des vectors singuliers 46 : C droites de MATA 47 : C 48 : C MATV = pointeur de la matrice (NVOIMI-1,NVOIMI-1) 49 : C des vectors singuliers gauches de MATA 50 : C 51 : C MATCEL = pointeur sur une matrice de travail 52 : C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient 53 : C momentanément b 54 : C la matrice (NTOTP-1,NVOIMI-1), transposé 55 : C de la pseudo-inverse de MATA 56 : C 57 : C MLRSIG = pointeur sur une liste de reels (NVOIMI-1) 58 : C qui contient lesvaleurs singuliers de MATA 59 : C 60 : C MLRCEL = pointeur sur une liste de reels 61 : C (NBMAXx(NBMAX-1)) (vecteur de travail) 62 : C 63 : C 64 : C SORTIES : COEF = pointeur sur une matrice 65 : C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient 66 : C la matrice (NTOTP,NVOIMI-1) pour calculer 67 : C le coefficients 68 : C 69 : C INDER = 1 -> probleme mal posé 70 : C 71 : C INDER = 2 -> erreur anormale 72 : C 73 : C----------------------------------------------------------------- 74 : C 75 : C HISTORIQUE (Anomalies et modifications eventuelles) 76 : C 77 : C HISTORIQUE : 10.05.2000, Creation 78 : C 79 : C----------------------------------------------------------------- 80 : C 81 : C 82 : IMPLICIT INTEGER(I-N) 83 : -INC CCOPTIO 84 : -INC SMCOORD 85 : -INC SMELEME 86 : -INC SMLENTI 87 : -INC SMLREEL 88 : C 89 : C----------------------------------------------------------------- 90 : C 91 : POINTEUR MLEVTO.MLENTI 92 : POINTEUR MLRSIG.MLREEL, MLRCEL.MLREEL 93 : C 94 : SEGMENT MATRIX 95 : REAL*8 MAT(N1,N2) 96 : ENDSEGMENT 97 : C 98 : POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX, MATCEL.MATRIX, 99 : & COEF.MATRIX 100 : C 101 : C**** Variables de COOPTIO 102 : C 103 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 104 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 105 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 106 : C & ,IECHO, IIMPI, IOSPI 107 : C & ,IDIM 108 : CC & ,MCOORD 109 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 110 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 111 : C & ,NORINC,NORVAL,NORIND,NORVAD 112 : C & ,NUCROU, IPSAUV 113 : C 114 : INTEGER NBMAX, NTOTP, NVOIMI, INDER, NOCENE, ICOORX, NOCEN1, IVOIS 115 : & , ILIG, ICOL, IVECT, IERSVD 116 : REAL*8 CONTOL, HELE, XCEN, YCEN, XCEN1, YCEN1, DX, DY, DX2, DY2 117 : & , HELE2, SVMAX, SVMIN 118 : C 119 : C 120 : INDER = 0 121 : C 122 : NOCENE = MLEVTO.LECT(1) 123 : HELE = 0.0D0 124 : ICOORX = ((IDIM + 1) * (NOCENE - 1))+1 125 : XCEN = MCOORD.XCOOR(ICOORX) 126 : YCEN = MCOORD.XCOOR(ICOORX+1) 127 : C 128 : C**** La submatrice de MATA qui nous interesse est (NTOTP-1)x(NVOIMI-1) 129 : C Amelioration: ca serait mieux travailler avec sa transposé 130 : C Mais invsvd.eso a ete teste avec m > n 131 : C 132 : DO IVOIS = 1, NTOTP-1, 1 133 : NOCEN1 = MLEVTO.LECT(IVOIS+1) 134 : ICOORX = ((IDIM + 1) * (NOCEN1 - 1))+1 135 : XCEN1 = MCOORD.XCOOR(ICOORX) 136 : YCEN1 = MCOORD.XCOOR(ICOORX+1) 137 : DX = XCEN1 - XCEN 138 : DY = YCEN1 - YCEN 139 : DX2 = DX * DX 140 : DY2 = DY * DY 141 : HELE = MAX(HELE,SQRT(DX2+DY2)) 142 : MATA.MAT(IVOIS,1) = DX 143 : MATA.MAT(IVOIS,2) = DY 144 : MATA.MAT(IVOIS,3) = 0.5D0 * DX2 145 : MATA.MAT(IVOIS,4) = 0.5D0 * DY2 146 : MATA.MAT(IVOIS,5) = DX * DY 147 : ENDDO 148 : C 149 : C******* On doit adimensioner MATA 150 : C 151 : HELE2 = HELE * HELE 152 : DO IVOIS = 1, NTOTP - 1, 1 153 : MATA.MAT(IVOIS,1) = MATA.MAT(IVOIS,1) / HELE 154 : MATA.MAT(IVOIS,2) = MATA.MAT(IVOIS,2) / HELE 155 : MATA.MAT(IVOIS,3) = MATA.MAT(IVOIS,3) / HELE2 156 : MATA.MAT(IVOIS,4) = MATA.MAT(IVOIS,4) / HELE2 157 : MATA.MAT(IVOIS,5) = MATA.MAT(IVOIS,5) / HELE2 158 : ENDDO 159 : C 159 : CALL INVSVD( NBMAX * (NBMAX -1),NTOTP-1,NVOIMI-1 ,MATA.MAT, 160 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD, 161 : & MLRCEL.PROG) 162 : C 163 : IF(IERSVD .NE. 0)THEN 164 : C OPEN(1,FILE='SVD.PBM') 165 : C WRITE(1,*) NBMAX * (NBMAX -1) 166 : C WRITE(1,*) NTOTP-1 167 : C WRITE(1,*) NVOIMI-1 168 : C WRITE(1,*) 169 : C DO IVOIS = 1, NTOTP - 1, 1 170 : C DO ILIG = 1, NVOIMI - 1 171 : C WRITE(1,*) MATA.MAT(IVOIS,ILIG) 172 : C ENDDO 173 : C ENDDO 174 : C CLOSE(1) 175 : INDER = 2 176 : GOTO 9999 177 : ENDIF 178 : C 179 : SVMAX = 0.0D0 180 : DO IVECT = 1, NVOIMI-1, 1 181 : SVMAX = MAX(SVMAX,MLRSIG.PROG(IVECT)) 182 : ENDDO 183 : C 184 : SVMIN = SVMAX 185 : DO IVECT = 1, NVOIMI-1, 1 186 : SVMIN = MIN(SVMIN,MLRSIG.PROG(IVECT)) 187 : ENDDO 188 : C 189 : IF((SVMIN / SVMAX) .LT. CONTOL)THEN 190 : DO IVOIS = 1 , NTOTP 191 : DO IVECT = 1, 5 192 : COEF.MAT(IVOIS,IVECT) = 0.0D0 193 : ENDDO 194 : ENDDO 195 : INDER = 1 196 : GOTO 9999 197 : ENDIF 198 : C 199 : C******* On doit calculer MATCEL (transpose de la pseudo inverse de MATA) 200 : C 201 : DO ICOL = 1, NVOIMI - 1, 1 202 : DO ILIG = 1, NTOTP - 1, 1 203 : MATCEL.MAT(ILIG,ICOL) = 0.0D0 204 : ENDDO 205 : ENDDO 206 : C 207 : DO IVECT = 1, NVOIMI - 1, 1 208 : DO ICOL = 1, NVOIMI - 1, 1 209 : DO ILIG = 1, NTOTP - 1, 1 210 : MATCEL.MAT(ILIG,ICOL) = MATCEL.MAT(ILIG,ICOL) + 211 : & (MATU.MAT(ILIG,IVECT) * MATV.MAT(ICOL,IVECT) 212 : & / MLRSIG.PROG(IVECT)) 213 : ENDDO 214 : ENDDO 215 : ENDDO 216 : C 217 : DO IVECT = 1, 5 218 : COEF.MAT(1,IVECT) = 0.0D0 219 : ENDDO 220 : C 221 : DO IVECT = 1, 2 222 : DO IVOIS = 1, NTOTP - 1, 1 223 : COEF.MAT(1+IVOIS,IVECT) = 224 : & MATCEL.MAT(IVOIS,IVECT) / HELE 225 : COEF.MAT(1,IVECT) = 226 : & COEF.MAT(1,IVECT) - 227 : & COEF.MAT(1+IVOIS,IVECT) 228 : ENDDO 229 : ENDDO 230 : C 231 : DO IVECT = 3, 5 232 : DO IVOIS = 1, NTOTP - 1, 1 233 : COEF.MAT(1+IVOIS,IVECT) = 234 : & MATCEL.MAT(IVOIS,IVECT) / HELE2 235 : COEF.MAT(1,IVECT) = 236 : & COEF.MAT(1,IVECT) - COEF.MAT(1+IVOIS,IVECT) 237 : ENDDO 238 : ENDDO 239 : 9999 RETURN 240 : C 241 : END 242 : 243 : 244 : 245 :
© Cast3M 2003 - All rights reserved.
Disclaimer