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