1 : C CGRESD SOURCE CHAT 05/01/12 21:53:31 5004 2 : SUBROUTINE CGRESD(KSTO, 3 : 1 NL,ILG, 4 : 2 IMAT,IPRC,IA,JA,KA, 5 : 3 IZB,IZP,ICOLD, 6 : 4 NPT,NPITE,NEFF,ICONV,EPI,IPOU,VPOU, 7 : 5 NIMPR,IPAT) 8 : C 9 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 10 : C C 11 : C RESOLUTION ITERATIVE D'UN SYSTEME C 12 : C SYMETRIQUE CREUX PAR LA METHODE DU C 13 : C GRADIENT CONJUGUE PRECONDITIONNE PAR LA C 14 : C DIAGONALE. C 15 : C C 16 : C***** INFOS SUR LA METHODE (LIGNE 0) C 17 : C C 18 : C DEUX MODES DE STOCKAGE SONT PROPOSES C 19 : C MORSE OU COMPRESSE. C 20 : C C 21 : C MORSE (KSTO=0) C 22 : C COMPRESSE (KSTO=1) C 23 : C C 24 : C POUR LE CAS DCG, ON RESOUD L'EQUATION SUIVANTE C 25 : C C 26 : C ( (K-1)*A*(K-1) )*xi = (K-1)*b C 27 : C C 28 : C AVEC K = D**(1/2) C 29 : C pi est la pseudo solution C 30 : C La solution x est obtenue en faisant C 31 : C x = (K-1)*xi C 32 : C C 33 : C***** INFOS SUR LES DIMENSIONS (LIGNE 1) C 34 : C C 35 : C NL (L'ORDRE DU SYSTEME) C 36 : C C 37 : C ILG (LONGUEUR LA POUR LE MORSE) C 38 : C (NOMBRE MAXI DE VALEURS PAR LIGNE DE LA C 39 : C MATRICE POUR LE COMPRESSE) C 40 : C C 41 : C***** INFOS SUR LA MATRICE (LIGNE 2) C 42 : C C 43 : C IMAT (LA MATRICE DANS UN LISTREEL) C 44 : C C 45 : C IPRC (LE PRECONDITIONNEMENT DANS UN LISTREEL) C 46 : C C 47 : C DCG : IPRC EST LE POINTEUR SUR UN LISTREEL C 48 : C CONTENANT D**(-1/2) CE LISTREEL A POUR C 49 : C DIMENSION NL AU MOINS C 50 : C C 51 : C IA (TABLEAU DE POINTEURS DU STOCKAGE MORSE) C 52 : C C 53 : C IL DOIT ETRE DIMENSIONNE A NL+1 AU MINIMUM C 54 : C C 55 : C JA (TABLEAU DE CONNECTIVITES DE LA MATRICE) C 56 : C C 57 : C IL EST STOCKE EN MORSE ET DE LONGUEUR LA C 58 : C C 59 : C KA (TABLEAU DE CONNECTIVITES EN MODE COMPRESSE) C 60 : C C 61 : C IL EST DE LA FORME KA(NL,NNZ) AVEC NNZ C 62 : C LE NOMBRE MAXI DE TERMES NON NULS POUR UNE C 63 : C LIGNE DE LA MATRICE C 64 : C C 65 : C NB : EN MORSE, KA N'EST PAS UTILISE C 66 : C EN COMPRESSE, IA ET JA NE SONT PAS UTILISES C 67 : C C 68 : C IA,JA ET KA SONT DES LISTENTI. C 69 : C C 70 : C***** INFOS SUR LE SECOND MEMBRE, LA SOLUTION INITIALE C 71 : C***** ET LA SOLUTION FINALE. (LIGNE 3) C 72 : C C 73 : C IZB SECOND MEMBRE EN ENTREE C 74 : C PAS MODIFIE. C 75 : C C 76 : C IZP ESTIMATION DE LA SOLUTION EN ENTREE C 77 : C LA SOLUTION EN SORTIE, IL FAUT LA C 78 : C CONSERVER POUR LA RESOLUTION SUIVANTE C 79 : C C 80 : C IZB ET IZP SONT DES LISTREEL DIMENSIONNES A NL C 81 : C C 82 : C ICOLD =1 (DEMARRAGE FROID) C 83 : C ICOLD!=1 (DEMARRAGE CHAUD A PARTIR DE LA C 84 : C SOLUTION INITIALE FOURNIE DANS IZP) C 85 : C C 86 : C***** INFOS SUR LA CONVERGENCE (LIGNE 4) C 87 : C C 88 : C NPT NOMBRE MAXI D'ITERATIONS EN ENTREE C 89 : C SI IL EST < 10, IL EST MIS A 10 C 90 : C NPITE FREQUENCE DES TESTS DE CONVERGENCE C 91 : C NEFF NOMBRE REEL EN SORTIE C 92 : C ICONV 0 NON CONVERGE, 1 CONVERGE C 93 : C EPI PRECISION DEMANDEE C 94 : C IPOU=1 LA PRECISION VA ETRE POUSSE D'UN FACTEUR C 95 : C VPOU. CECI PEUT S'AVERER UTILE LORS D'UN C 96 : C CALCUL ITERATIF OU ON VEUT DES LE DEBUT C 97 : C UNE SOLUTION PRECISE POUR EVITER DE PARTIR C 98 : C SUR UNE MAUVAISE PISTE. VPOU A UNE VALEUR C 99 : C PAR DEFAUT DE 100. C 100 : C C 101 : C***** INFOS SUR LE CONTROLE DES IMPRESSIONS (LIGNE 5) C 102 : C C 103 : C NIMPR =1 (INFOS SUR LA CONVERGENCE) C 104 : C NIMPR!=1 (SILENCIEUX SAUF EN CAS DE DIVERGENCE) C 105 : C C 106 : C IPAT = (FACULTATIF, PAR EXEMPLE UN PAS DE TEMPS)C 107 : C C 108 : C RSETD(B,A) B <- A C 109 : C ADIVEC(A,B,C,alpha) C <- A + alpha*B C 110 : C PIMDV(A,B,C) C <- A*B C 111 : C C 112 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 113 : C 114 : IMPLICIT INTEGER(I-N) 115 : IMPLICIT REAL*8 (A-H,O-Z) 116 : -INC SMLENTI 117 : POINTEUR IA.MLENTI 118 : POINTEUR JA.MLENTI 119 : POINTEUR KA.MLENTI 120 : -INC SMLREEL 121 : POINTEUR IMAT.MLREEL 122 : POINTEUR IPRC.MLREEL 123 : POINTEUR IZB.MLREEL 124 : POINTEUR IZBB.MLREEL 125 : POINTEUR IZP.MLREEL 126 : POINTEUR IZPS.MLREEL 127 : POINTEUR IZZPR.MLREEL 128 : POINTEUR IZPMV.MLREEL 129 : POINTEUR IZTMP.MLREEL 130 : POINTEUR IZDIR.MLREEL 131 : 132 : DATA VPOU0/100./ 133 : 134 : IF(KSTO.EQ.0) LA=ILG 135 : IF(KSTO.EQ.1) NNZ=ILG 136 : IF(IPOU.EQ.1) THEN 137 : VPOUE=VPOU0 138 : IF(VPOU.NE.VPOU0) VPOUE=VPOU 139 : ENDIF 140 : 140 : CALL OOOMRU(1) 141 : 142 : C On sauve le second membre en le dupliquant dans IZBB 143 : C qui lui sera ecrase. 144 : SEGACT IZB 145 : JG=IZB.PROG(/1) 146 : SEGINI IZBB 146 : CALL RSETD(IZBB.PROG,IZB.PROG,JG) 147 : SEGDES IZB 148 : 149 : SEGACT IZP 150 : 151 : JG=NL 152 : SEGINI IZPS 153 : SEGINI IZZPR 154 : SEGINI IZPMV 155 : SEGINI IZTMP 156 : SEGINI IZDIR 157 : 158 : SEGACT IMAT 159 : SEGACT IPRC 160 : IF(KSTO.EQ.0) SEGACT IA 161 : IF(KSTO.EQ.0) SEGACT JA 162 : IF(KSTO.EQ.1) SEGACT KA 163 : 164 : IF(NPT.LT.10) NPT=10 165 : IF(NPT.GT.NL) NPT=NL 166 : 167 : C Mettre le second membre et l'estimation initiale 168 : C a l'echelle, B <- D-(1/2)*B, PIM <- PIM/(D**(-1/2)) 168 : CALL PIMDV(IZBB.PROG,IPRC.PROG,IZBB.PROG,NL) 168 : CALL DIVISE(NL,IZP.PROG,IZP.PROG,IPRC.PROG) 168 : IF(ICOLD.EQ.1) CALL INITD (IZP.PROG,NL,0.D0) 168 : CALL RSETD(IZPS.PROG,IZP.PROG,NL) 168 : CALL PIMDV(IZP.PROG,IPRC.PROG,IZTMP.PROG,NL) 169 : IF(KSTO.EQ.0) THEN 169 : C CALL PMVM(A,IA,JA,X,Y,NL,LA) 169 : CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT, 170 : 1 IZTMP.PROG,IZPMV.PROG,NL,LA) 171 : ELSE 171 : C CALL PMVC(A,KA,X,Y,NL,NNZ) 171 : CALL PMVC(IMAT.PROG,KA.LECT, 172 : 1 IZTMP.PROG,IZPMV.PROG,NL,NNZ) 173 : ENDIF 173 : CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL) 173 : CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,-1.d0,NL) 174 : ANORB=SQRT(PISCAL(IZBB.PROG,IZBB.PROG,NL)) 174 : CALL RSETD(IZZPR.PROG,IZBB.PROG,NL) 174 : CALL RSETD(IZDIR.PROG,IZZPR.PROG,NL) 175 : SCAR1=PISCAL(IZBB.PROG,IZDIR.PROG,NL) 176 : C 177 : C===========DEBUT DE LA BOUCLE ITERATIVE=====600======== 178 : C 179 : NEFF=0 180 : NCOMP=0 181 : ICONV=0 182 : DO 600 I=1,NPT 183 : NEFF=NEFF+1 184 : NCOMP=NCOMP+1 184 : CALL PIMDV(IZDIR.PROG,IPRC.PROG,IZTMP.PROG,NL) 185 : IF(KSTO.EQ.0) THEN 185 : C CALL PMVM(A,IA,JA,X,Y,NL,LA) 185 : CALL PMVM(IMAT.PROG,IA.LECT,JA.LECT, 186 : 1 IZTMP.PROG,IZPMV.PROG,NL,LA) 187 : ELSE 187 : C CALL PMVC(A,KA,X,Y,NL,NNZ) 187 : CALL PMVC(IMAT.PROG,KA.LECT, 188 : 1 IZTMP.PROG,IZPMV.PROG,NL,NNZ) 189 : ENDIF 189 : CALL PIMDV(IZPMV.PROG,IPRC.PROG,IZPMV.PROG,NL) 190 : SCADIR=PISCAL(IZDIR.PROG,IZPMV.PROG,NL) 191 : ALPHA=SCAR1/SCADIR 191 : CALL ADIVEC(IZP.PROG,IZDIR.PROG,IZP.PROG,ALPHA,NL) 192 : ALP=-1.D0*ALPHA 192 : CALL ADIVEC (IZBB.PROG,IZPMV.PROG,IZBB.PROG,ALP,NL) 192 : CALL RSETD(IZZPR.PROG,IZBB.PROG,NL) 193 : SCAR2=PISCAL(IZBB.PROG,IZZPR.PROG,NL) 194 : BETA=SCAR2/SCAR1 194 : CALL ADIVEC (IZZPR.PROG,IZDIR.PROG,IZDIR.PROG,BETA,NL) 195 : IF (NCOMP.EQ.NPITE.OR.I.EQ.NPT) THEN 196 : INDIC=0 197 : IF (I.EQ.NPT) INDIC=1 198 : C On pousse la precision si necessaire 199 : ZEPS=EPI 200 : IF(IPOU.EQ.1) ZEPS=EPI/VPOUE 200 : CALL COTEST (IZP.PROG,IZPS.PROG, 201 : 1 ZEPS,DELTA,DELTAP, 202 : 2 INDIC,ICONV,NL) 203 : SEPA=DELTA 204 : RSEPA=DELTA/DELTAP 204 : CALL RSETD(IZPS.PROG,IZP.PROG,NL) 205 : NCOMP=0 206 : ENDIF 207 : IF (ICONV.EQ.1) GOTO 601 208 : SCAR1=SCAR2 209 : 600 CONTINUE 210 : 601 CONTINUE 210 : CALL MTABL(IZBB.PROG,RES,IMA,NL) 211 : RRES=RES/ANORB 212 : C Pour le DCG la solution est D**(-1/2)*IZP 213 : C Pour le DCG on remet dans PIM, la solution 213 : CALL PIMDV(IZP.PROG,IPRC.PROG,IZP.PROG,NL) 214 : IF (ICONV.NE.1) THEN 215 : WRITE(6,1111)IPAT,NPT,RSEPA,DELTAP,RES,IZBB.PROG(IMA) 216 : ELSE 217 : IF(NIMPR.EQ.1) WRITE(6,2222)IPAT,EPI, 218 : 1 RES,SEPA,I, 219 : 2 RRES,RSEPA,IMA 220 : ENDIF 221 : 222 : SEGDES IZP 223 : 224 : SEGSUP IZBB 225 : SEGSUP IZPS 226 : SEGSUP IZZPR 227 : SEGSUP IZDIR 228 : SEGSUP IZTMP 229 : SEGSUP IZPMV 230 : 231 : SEGDES IMAT*(NOMOD,MRU) 232 : SEGDES IPRC*(NOMOD,MRU) 233 : IF(KSTO.EQ.0) SEGDES IA*(NOMOD,MRU) 234 : IF(KSTO.EQ.0) SEGDES JA*(NOMOD,MRU) 235 : IF(KSTO.EQ.1) SEGDES KA*(NOMOD,MRU) 236 : 236 : CALL OOOMRU(0) 237 : 238 : C 239 : 1111 FORMAT(1X,'*CGRESD* ',I6,' NON CONV EN ',I5, 240 : 1' ITERA, RSEPA = ',D8.2,1X,'DELTAP = ',D8.2,1X, 241 : 2' RES = ',D8.2,' P = ',D8.2) 242 : 2222 FORMAT(1X,'*CGRESD* ',I6,' *PREC ',D8.2, 243 : 1 ' * RESIDU-SEPARA @ NITE ',D8.2,1X,D8.2,1X,I5, 244 : 2 ' * RRES-RSEP ',D8.2,1X,D8.2, 245 : 3 ' * ELEMENT ',I7) 246 : C 247 : C FIN RESOLUTION 248 : C 249 : RETURN 250 : END 251 : 252 :
© Cast3M 2003 - All rights reserved.
Disclaimer