C PRBCG2 SOURCE GOUNAND 11/07/21 21:15:31 7046 SUBROUTINE PRBCG2(KMORS,KISA,KS2B,MATRIK,MAPREC,LRES,LNMV, $ INCX,ITER,INMV, $ RESID,KPREC, $ BRTOL,ICALRS,IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : PRBCG2 C DESCRIPTION : C Préparation de la résolution d'un système linéaire Ax=b C par une méthode de bi-gradient conjugué stabilisé C (Bi-CGSTAB <~> Bi-CG + GMRES(1)) C préconditionnée ou non. C C A peut etre symétrique (mais mieux vaut utiliser C le gradient conjugué) C ou non symétrique (mais la méthode peut ne pas converger C voir paramètre BRTOL en particulier). C C Les 4 préconditionnements disponibles pour cette méthode C sont : C * Jacobi (diagonal) C * D-ILU (Diagonal Incomplete LU factorization) C * ILU(0) (Incomplete LU factorization C of level zero ie Choleski incomplet) C * MILU(0) (Modified ILU(0)) relaxé C 5 : préconditionnement ILUT (dual truncation) C 6 : préconditionnement ILUT2 (une variante du C précédent qui remplit mieux la mémoire et C fonctionne mieux quelquefois) C C Ces préconditionneurs sont supposés déjà construits (ie C factorisés), voir les subroutines MEIDIA, PRDILU, PRILU0, C PRMILU et PRILUT. C C Ce sous-programme est en fait une interface aux sous-programmes : C bcgs, bcgsd, bcgsdi, bcgsi0 C qui sont en Fortran presque pur (pour raison de rapidité) C et effectuent la résolution du système linéaire C proprement dite. C C LANGAGE : ESOPE C AUTEUR : Stéphane GOUNAND (CEA/DRN/DMT/SEMT/TTMF) C mél : gounand@semt2.smts.cea.fr C REFERENCE (bibtex-like) : C @BOOK{templates, C AUTHOR={R.Barrett, M.Berry, T.F.Chan, J.Demmel, J.Donato, C J.Dongarra, V.Eijkhout, R.Pozo, C.Romine, C H. Van der Vorst}, C TITLE={Templates for the Solution of Linear Systems : C Building Blocks for Iterative Methods}, C PUBLISHER={SIAM}, YEAR={1994}, ADDRESS={Philadelphia,PA} } C -> URL : http://www.netlib.org/templates/Templates.html C*********************************************************************** C APPELES (Calcul) : BCGS, BCGSD, BCGSDI, BCGSI0 C*********************************************************************** C ENTREES : MATRIK, MAPREC, KPREC, BRTOL, IMPR C ENTREES/SORTIES : INCX, ITER, RESID C SORTIES : LRES, IRET C CODE RETOUR (IRET) : 0 si ok C <0 si problème C >0 si l'algorithme appelé n'a pas convergé C MATRIK : pointeur sur segment MATRIK de l'include SMMATRIK C on pioche dedans les informations nécessaires C (différents pointeurs, nb. de ddl...) C MAPREC : pointeur sur segment MATRIK de l'include SMMATRIK C on va l'utiliser comme préconditionneur. C KPREC : type de préconditionnement à effectuer C 0 : pas de préconditionnement C 1 : préconditionnement Jacobi (diagonal) C 2 : préconditionnement D-ILU C 3 : préconditionnement ILU(0) (Choleski incomplet) C 4 : préconditionnement MILU(0) C (Choleski incomplet modifié) C 5 : préconditionnement ILUT (dual truncation) C 6 : préconditionnement ILUT2 (une variante du C précédent qui remplit mieux la mémoire et C fonctionne mieux quelquefois) C BRTOL : REAL*8 C 'Breakdown tolerance' pour les méthodes de type C BiCGSTAB. Si un certain produit scalaire de vecteurs C "direction" est inférieur à cette tolérance, la C méthode ne peut pas continuer. C IMPR : niveau d'impression (0..4) C INCX : pointeur sur segment IZA de l'include SMMATRIK C C'est le vecteur des inconnues. C En entrée, contient l'estimation x(0). C En sortie, la solution trouvée (si la méthode C a convergé). C ITER : type INTEGER. C En entrée, contient le nombre maximum C d'itérations à effectuer. C En sortie, contient le nombre d'itérations C réellement effectuées. C RESID : type REAL*8. C En entrée, la valeur maximale du résidu à C convergence de l'algorithme mesurée comme suit : C norme[L2](b - A*x) / norme[L2]( b ) C En sortie, la valeur effective de cette mesure. C LRES : pointeur sur segment MLREEL (include SMLREEL) C contient l'historique de la valeur de RESID en C fonction du nombre d'itérations effectuées. C C*********************************************************************** C VERSION : v1, 02/04/98, version initiale C HISTORIQUE : v1, 02/04/98, création C HISTORIQUE : 09/02/99 : ajout préconditionneur <> MATRIK C HISTORIQUE : 20/12/99 : interfaçage avec le nouveau bcgsi0 C (factorisations incomplètes stockées au format MSR (Modified Sparse C Row) (voir la doc de Sparskit version 2+) C http://www.cs.umn.edu/Research/arpa/SPARSKIT/sparskit.html C HISTORIQUE : 22/03/00 : ajout des préconditionneurs ILUT C HISTORIQUE : 09/04/04 : ajout des préconditionneurs ILUTP C HISTORIQUE : 10/02/06 : Nettoyage C HISTORIQUE : C*********************************************************************** C Prière de PRENDRE LE TEMPS de compléter les commentaires C en cas de modification de ce sous-programme afin de faciliter C la maintenance ! C*********************************************************************** * * .. Includes et pointeurs associés .. -INC CCOPTIO -INC SMLREEL INTEGER JG POINTEUR LRES.MLREEL -INC SMLENTI POINTEUR LNMV.MLENTI -INC SMMATRIK POINTEUR MAPREC.MATRIK POINTEUR KMORS.PMORS POINTEUR KISA.IZA POINTEUR KS2B.IZA POINTEUR INCX.IZA POINTEUR INVDIA.IZA POINTEUR ILUM.PMORS POINTEUR ILUI.IZA * .. Parameters * This is a breakdown tolerance for BiCGSTAB-type method REAL*8 BRTOL * .. Work Vectors for BCGS.. POINTEUR IR.IZA,IRTLD.IZA,IP.IZA,IPHAT.IZA POINTEUR IV.IZA,IS.IZA,ISHAT.IZA,IT.IZA *STAT-INC SMSTAT * .. Scalar Arguments .. INTEGER ITER, KPREC, IMPR, IRET REAL*8 RESID * .. * Vars reqd for STOPTEST2 * REAL*8 TOL, BNRM2 * .. * .. External subroutines .. * EXTERNAL STOPTEST2 INTEGER NBVA,NJA,NTT,NTTPRE * .. * .. Executable Statements .. * IRET = 0 * * On récupère les paramètres * SEGACT MATRIK SEGACT MAPREC IF (KSYM.EQ.0) THEN IF (IMPR.GT.0) THEN WRITE(IOIMP,*) 'MATRIK',MATRIK,'symétrique : ', $ 'Use a Conjugate Gradient instead !' C IRET=-2 C GOTO 9999 ENDIF ENDIF * Pour le préconditionneur ILUM =MAPREC.KIDMAT(6) ILUI =MAPREC.KIDMAT(7) IDMAT=MAPREC.KIDMAT(1) SEGACT IDMAT INVDIA=IDIAG SEGDES IDMAT SEGACT KMORS NTT =KMORS.IA(/1)-1 * NJA =KMORS.JA(/1) SEGACT KISA SEGACT KS2B SEGACT INCX*MOD IF (KPREC.NE.0) THEN C Paramètres des préconditionnements diagonaux et D-ILU IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN C Est-il compatible avec la matrice ? SEGACT INVDIA NTTPRE=INVDIA.A(/1) IF (NTTPRE.NE.NTT) THEN WRITE(IOIMP,*) 'La matrice et son préconditionnement' WRITE(IOIMP,*) 'ne sont pas compatibles...' WRITE(IOIMP,*) 'NTT=',NTT WRITE(IOIMP,*) 'NTTPRE=',NTTPRE IRET=-2 GOTO 9999 ENDIF C Paramètres des préconditionnements ILU(0), MILU(0), ILUT, ILUT2 C ,ilutp, ilutpg, ilutpg2 ELSEIF (KPREC.GE.3.AND.KPREC.LE.10) THEN SEGACT ILUM SEGACT ILUI * SEGPRT,ILUM * SEGPRT,ILUI NTTPRE=ILUM.IA(/1) IF (NTTPRE.NE.NTT) THEN WRITE(IOIMP,*) 'La matrice et son préconditionnement', $ 'ne sont pas compatibles...' WRITE(IOIMP,*) 'NTT=',NTT,' NTTPRE=',NTTPRE IRET=-2 GOTO 9999 ENDIF ELSE INVDIA=0 WRITE(IOIMP,*) 'Erreur de programmation' GOTO 9999 ENDIF ENDIF C C Initialisation de l'historique de convergence C JG=ITER+1 SEGINI LNMV SEGINI LRES C C C Initialisation des vecteurs de travail pour BiCGSTAB C *STAT CALL INMSTA(MSTAT,0) NBVA=NTT SEGINI IR,IRTLD,IP,IPHAT SEGINI IV,IS,ISHAT,IT *STAT CALL PRMSTA(' PRBCG2 : objets temporaires',MSTAT,1) *STAT CALL SUMSTA(MSTAT,0) C C Appel du BiCGSTAB C CALL BCG2(KMORS,KISA,KS2B,INCX, $ KPREC,INVDIA,ILUM,ILUI, $ LRES,LNMV, $ IR,IRTLD,IP,IPHAT, $ IV,IS,ISHAT,IT, $ ITER,INMV,BRTOL,RESID,ICALRS,IMPR,IRET) C C Désactivation C SEGSUP IR,IRTLD,IP,IPHAT SEGSUP IV,IS,ISHAT,IT JG=ITER+1 SEGADJ LRES SEGDES LRES SEGADJ LNMV SEGDES LNMV IF (KPREC.EQ.1.OR.KPREC.EQ.2) THEN SEGDES INVDIA ELSEIF (KPREC.GE.3.AND.KPREC.LE.9) THEN SEGDES ILUI SEGDES ILUM ENDIF SEGDES INCX SEGDES KS2B SEGDES KISA SEGDES KMORS SEGDES MAPREC SEGDES MATRIK C C A breakdown has occured in the BiCGSTAB method C IF (IRET.LT.0) GOTO 9999 * * Normal termination * RETURN * * Format handling * * 1002 FORMAT(10(1X,1PE11.4)) * * Error handling * 9999 CONTINUE WRITE(IOIMP,*) 'An error was detected in prbcg2.eso' RETURN * * End of PRBCG2 * END