C PRAGMG SOURCE GOUNAND 25/04/30 21:15:28 12258 SUBROUTINE PRAGMG(KMORS,KIZA,KS2B,IPBLOC,KTYPI,LRES, $ LNMV,INCX,ITER,INMV, $ RESID,KPREC, $ NRESTS,ICALRS,IDDOT,IMVEC,KTIME,LTIME,LDUMP,ISMOOT, $ IMPR,IRET) IMPLICIT REAL*8 (A-H,O-Z) IMPLICIT INTEGER (I-N) C*********************************************************************** C NOM : PRAGMG C DESCRIPTION : C Préparation de la résolution d'un système linéaire Ax=b C par une méthode Multigrille Algébrique 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*********************************************************************** C*********************************************************************** C VERSION : v1, 17/06/08, version initiale C HISTORIQUE : 17/06/08 : Creation 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 PPARAM -INC CCOPTIO -INC CCREEL -INC SMLREEL INTEGER JG POINTEUR LRES.MLREEL -INC SMLENTI POINTEUR LNMV.MLENTI POINTEUR IBLOCK.MLENTI POINTEUR IPBLOC.MLENTI POINTEUR IPERM.MLENTI POINTEUR JPERM.MLENTI -INC SMMATRIK POINTEUR KMORS.PMORS POINTEUR KIZA.IZA,AIZA.IZA POINTEUR KS2B.IZA,RES.IZA POINTEUR INCX.IZA,INCDX.IZA -INC SMTABLE POINTEUR KTIME.MTABLE CHARACTER*8 CHARI LOGICAL LTIME,LOGII,LBLOCK,LOK,LDUMP,LSIGN * .. Parameters * This is a breakdown tolerance for BiCGSTAB-type method REAL*8 BRTOL *STAT-INC SMSTAT * .. Scalar Arguments .. INTEGER ITER, KPREC, IMPR, IRET INTEGER ICNT REAL*8 RESID,TOL,BNRM2,RNRM2,XFACT * .. * Vars reqd for STOPTEST2 * REAL*8 TOL, BNRM2 * .. * .. External subroutines .. * EXTERNAL STOPTEST2 INTEGER NBVA,NJA,NTT,NTTPRE REAL*8 GNRM2 EXTERNAL GNRM2 * .. * .. Executable Statements .. * * WRITE(IOIMP,*) 'Debut de pragmg' IRET = 0 * * On récupère les paramètres * * segact kmors segact kiza NTT =KMORS.IA(/1)-1 NJA =KMORS.JA(/1) NBVA=KIZA.A(/1) SEGINI,RES=KS2B nbva=INCX.A(/1) SEGINI INCDX IBLOCK=IPBLOC NBLOCK=IBLOCK.LECT(/1)-1 C C Initialisation de l'historique de convergence C JG=ITER+1 SEGINI LNMV SEGINI LRES * * Autres paramètres * IJOB=0 IF (LDUMP) THEN IPRINT=-978 ELSE IF (IMPR.GT.2) THEN IPRINT=IOIMP ELSE IPRINT=-1 ENDIF ENDIF TOL=RESID ICNT=0 * * AGMG ne fait que du résidu relatif d'où les petits ajustements * suivants. * * res = b - AX CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,RES) BNRM2=GNRM2(KS2B) RNRM2=GNRM2(RES) WRITE(IOIMP,*) 'RNRM2=',RNRM2 IF (ICALRS.EQ.1) BNRM2=RNRM2 IF (BNRM2.LT.XPETIT) BNRM2=1.D0 RESID=RNRM2 / BNRM2 ICNT=ICNT+1 LRES.PROG(ICNT)=RESID LNMV.LECT(ICNT)=1 * IF (RESID.LE.TOL) THEN ITER=0 GOTO 30 ENDIF * IF (ICALRS.EQ.0) THEN XFACT=1.D0/RESID TOL=TOL*XFACT ENDIF * * Changement automatique du signe des lignes de la matrice * et du second membre si le terme diagonal est négatif. * SEGACT,KMORS segact,KIZA * LSIGN=.FALSE. DO I=1,NTT DO J=KMORS.IA(I),(KMORS.IA(I+1)-1) IF (KMORS.JA(J).EQ.I) THEN IF (KIZA.A(J).LT.XZERO) THEN LSIGN=.TRUE. IDEB=I GOTO 20 ELSE LSIGN=.FALSE. ENDIF ENDIF ENDDO ENDDO 20 CONTINUE * IF (LSIGN) THEN NSIGN=0 SEGINI,AIZA=KIZA DO I=IDEB,NTT IFOUND=0 DO J=KMORS.IA(I),(KMORS.IA(I+1)-1) * WRITE(IOIMP,*) 'I,J=',I,J * WRITE(IOIMP,*) 'JA(J)=',KMORS.JA(J) IF (KMORS.JA(J).EQ.I) THEN * WRITE(IOIMP,*) 'AISA.A(J)=',AISA.A(J) IF (AIZA.A(J).LT.XZERO) THEN IFOUND=-1 GOTO 10 ENDIF ENDIF ENDDO 10 CONTINUE IF (IFOUND.EQ.-1) THEN NSIGN=NSIGN+1 RES.A(I)=-1.D0*RES.A(I) DO J=KMORS.IA(I),(KMORS.IA(I+1)-1) AIZA.A(J)=-1.D0*AIZA.A(J) ENDDO ENDIF ENDDO WRITE(IOIMP,*) NSIGN,' matrix rows changed sign.' ELSE AIZA=KIZA ENDIF * * To use the standard AGMG library, you should make sure that it * is compiled with same compiler and options than Castem's and that * the sequential example furnished with AGMG works. * Then it is sufficient to put agmg_seq.o and agmg.o in the * current directory and use essai for linking. * * No interface to the parallel version of agmg for now * WRITE(IOIMP,*) '***********************************' WRITE(IOIMP,*) ' ' WRITE(IOIMP,*) 'Please uncomment the CALL AGMG(...', $ ' in the pragmg.eso subroutine if you wish to use' WRITE(IOIMP,*) 'the AGMG library by Y. Notay' WRITE(IOIMP,*) ' ' WRITE(IOIMP,*) '***********************************' * 251 2 *Tentative d'utilisation d'une option non implémentée CALL ERREUR(251) IRET=-1 GOTO 9999 ********************************************************************* * * AGMG v3.3.7_block call (standard sequential version) * * WRITE(IOIMP,*) '4' * CALL ECMORS(KMORS,AIZA,3) c$$$ c$$$ ITR2=ITER-1 c$$$ instance=0 c$$$ NREST=NRESTS c$$$ IF (KTYPI.EQ.7.OR.KTYPI.EQ.8) then c$$$ call DAGMG_INTPARAM('smoothtype',ISMOOT,instance) c$$$ CALL DAGMG(NTT,AIZA.A,KMORS.JA,KMORS.IA,RES.A,INCDX.A, c$$$ $ IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT) c$$$ CALL dagmg_status(istat,instance) c$$$ CALL dagmg_niter(itr2,instance) c$$$ ELSE c$$$ if(ktypi.eq.10) then c$$$ call DAGMG_STOKES_INTPARAM('smoothtype',-1,instance) c$$$ endif c$$$ if(ktypi.eq.11) then c$$$ call DAGMG_STOKES_INTPARAM('smoothtype',-51,instance) c$$$ call DAGMG_STOKES_INTPARAM('pointcoarsening',1,instance) c$$$ call DAGMG_STOKES_INTPARAM('Ltransform',1,instance) c$$$ endif c$$$ CALL DAGMG_STOKES(NTT,AIZA.A,KMORS.JA,KMORS.IA,RES.A,INCDX.A, c$$$ $ IJOB,IPRINT,NREST,ITR2,TOL,NBLOCK,IBLOCK.LECT) c$$$ CALL dagmg_stokes_status(istat,instance) c$$$ CALL dagmg_stokes_niter(itr2,instance) c$$$ ENDIF c$$$* X = X_0 + \delta X c$$$ CALL GAXPY(1.D0,INCDX,INCX) c$$$* ON peut l'eviter car deja calcule par AGMG c$$$** Résidu final c$$$* CALL GCOPY(KS2B,RES) c$$$* CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,RES) c$$$* RNRM2=GNRM2(RES) c$$$* RESID=RNRM2/BNRM2 c$$$* c$$$ IF (istat.NE.0) THEN c$$$ IRET=1 c$$$ ELSE c$$$ IRET=0 c$$$ ENDIF c$$$ ITER=ITR2 c$$$* c$$$* write(ioimp,*) 'ITER,nbva=',ITER,nbva c$$$ JMAX=MIN(ITER+4,nbva) c$$$ DO J=5,JMAX c$$$ ICNT=ICNT+1 c$$$* write(ioimp,*) 'ICNT=',ICNT c$$$ RESID=RES.A(J)/BNRM2 c$$$ LRES.PROG(ICNT)=RESID c$$$ LNMV.LECT(ICNT)=LNMV.LECT(ICNT-1)+1 c$$$ ENDDO c$$$ INMV=1+ITER * * End of AGMG call (standard sequential version) * ********************************************************************* * IF (LTIME) THEN * CHARI='MGAGGREG' * CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, * $ 'ENTIER ',IAT,XVALR,CHARR,LOGIR,IRETR) * CHARI='MGSOLUTI' * CALL ECCTAB(KTIME,'MOT ',IVALI,XVALI,CHARI,LOGII,IRETI, * $ 'ENTIER ',IST,XVALR,CHARR,LOGIR,IRETR) ENDIF * C C Désactivation C 30 CONTINUE * segprt,lres * segprt,lnmv JG=ICNT SEGADJ LRES SEGDES LRES SEGADJ LNMV SEGDES LNMV * segprt,lres * segprt,lnmv SEGSUP INCDX SEGDES INCX SEGDES KS2B SEGSUP RES IF (LSIGN) SEGSUP AIZA SEGDES KIZA SEGDES KMORS C C A breakdown has occured in the CGS 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 pragmg.eso' RETURN * * End of PRAGMG * END