C GRQUA4 SOURCE CHAT 05/01/13 00:21:35 5004 SUBROUTINE GRQUA4(LISTPO,MLEVOI,MLEVO2,MLRDIS, & KRICEN,CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO, & MATA,MATU,MATV,MATCEL,MLRSIG,MLRCEL,COEF,INDER) C------------------------------------------------------------------ C PROJET : CASTEM 2000 C C NOM : GRQUA4 C C DESCRIPTION : Augmente le stencil pour le calcul du gradient. C C LANGAGE : FORTRAN 77 + ESOPE 2000 C C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/LTMF C C----------------------------------------------------------------- C C APPELES (E/S) : C C APPELES (Calcul): C C----------------------------------------------------------------- C C ENTREES : LISTPO = pointeur sur une liste de point contenant C les centres et ses voisins (voir la C description du segment) C C MLEVOI = pointeur sur une liste MLENTI C (dimensions = nombre total de points) C C MLEVO2 = pointeur sur une liste MLENTI qui va contenir C les points "voisins des voisins" de centre C C MLRDIS = pointeur sur une liste MLREEL qui va contenir C les distances des points MLEVO2 du point C centre C C KRICEN = pointeur sur une liste MLENTI qui contient C la correpondance LOCAL/GLOBAL du maillage C des points centres C C CONTOL = REAL 8 t.q. 1/cond(MATA) > CONTOL C Si non rank(MATA) < NVOIMI - 1 C C NBMAX, NTOTP, NVOIMI : INTEGERs qui interviennent C dans les dimensions de C certains matrices C N.B. NTOTP = E/S C C C MATA = pointeur d'une matrice C (NBMAXx(NBMAX-1),NVOIMI-1) qui contient la C matrice (NTOTP-1,NVOIMI-1), a pseudoinverser C C MATU = pointeur d'une matrice C (NBMAXx(NBMAX-1),NVOIMI-1) on met la matrice C (NTOTP-1,NVOIMI-1) des vectors singuliers C droites de MATA C C MATV = pointeur de la matrice (NVOIMI-1,NVOIMI-1) C des vectors singuliers gauches de MATA C C MATCEL = pointeur sur une matrice de travail C (NBMAX (NBMAX - 1),NVOIMI-1) qui contient C momentanément b C la matrice (NTOTP-1,NVOIMI-1), transposé C de la pseudo-inverse de MATA C C MLRSIG = pointeur sur une liste de reels (NVOIMI-1) C qui contient lesvaleurs singuliers de MATA C C MLRCEL = pointeur sur une liste de reels C (NBMAXx(NBMAX-1)) (vecteur de travail) C C E/S : NTOTP = nombre total des voisins du point centre C dont on veut augmenter le stencil C C MLEVTO = pointeur sur un list d'entiers C (NBMAX - 1) * NBMAX : C MLEVTO.LECT(1) = point ou je veux calculer C gradient et hessien C MLEVTO.LECT(i), i=2, NTOTP = voisins C C SORTIES : COEF = pointeur sur la matrice qui contient le C coefficient pour le calcul de l'hessien C C INDER = 1 -> probleme mal posé C = 2 -> erreur anormale C C----------------------------------------------------------------- C C HISTORIQUE (Anomalies et modifications eventuelles) C C HISTORIQUE : 11.05.2000, Creation C C----------------------------------------------------------------- C C IMPLICIT INTEGER(I-N) -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMLENTI -INC SMLREEL C C----------------------------------------------------------------- C POINTEUR MLEVOI.MLENTI, MLEVO2.MLENTI, MLRDIS.MLREEL, & KRICEN.MLENTI, MLEVTO.MLENTI, MLRSIG.MLREEL, & MLRCEL.MLREEL C SEGMENT MATRIX REAL*8 MAT(N1,N2) ENDSEGMENT C POINTEUR MATA.MATRIX, MATU.MATRIX, MATV.MATRIX, MATCEL.MATRIX, & COEF.MATRIX C SEGMENT LISTPO INTEGER LISPO1(NPOILI),INDEX1(NBCEN+1) ENDSEGMENT C C**** LISTPO: list de points centres avec leur voisins C et position: C C LISPO1(k) t.q. C LISPO1(INDEX1(j)) = numero global du centre j C LISPO1(INDEX1(j)),...,LISPO1(INDEX1(j+1))-1 = voisins du centre C LISPO1(INDEX1(j)) C C C**** Variables de COOPTIO C C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES C & ,IECHO, IIMPI, IOSPI C & ,IDIM CC & ,MCOORD C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU C & ,NORINC,NORVAL,NORIND,NORVAD C & ,NUCROU, IPSAUV C INTEGER NVOIMI, NTOTP, INDER, INEUVO, NOCENE, LAST0, IVOIS, NOCEN1 & , ICEN, IPOS1, NVOIS1, NOCEN2, ICOORX, IVOIS2, IVOIS3, IPOS0 & , NBMAX C REAL*8 XCEN, YCEN, ZCEN, DIST, DIST1, CONTOL C C INEUVO = nombre de nouveaux voisins C INEUVO = 0 C INDER = 0 NOCENE = MLEVTO.LECT(1) C ICOORX = ((IDIM + 1) * (NOCENE - 1))+1 XCEN = MCOORD.XCOOR(ICOORX) YCEN = MCOORD.XCOOR(ICOORX+1) IF(IDIM .EQ. 3) ZCEN = MCOORD.XCOOR(ICOORX+2) C LAST0 = NOCENE MLEVOI.LECT(NOCENE) = -1 C C**** On mets -1 dans le dernier elt de la liste chaîné C -1 != 0 -> on sait que NOCENE fait partir de l'elt de C noeuds + voisins qu'on veut realizer! C DO IVOIS = 2, NTOTP, 1 NOCEN1 = MLEVTO.LECT(IVOIS) MLEVOI.LECT(NOCEN1) = LAST0 LAST0 = NOCEN1 ENDDO C C**** Liste chainé du noeud et de ses principaux voisins remplie! C NOCEN1 = LAST0 C C**** Si NOCEN1 = NOCENE, ça veut dire que le pouvre est tout C seul! C'est tres triste! C IF(NOCEN1 .EQ. NOCENE) THEN MOTERR(1:40) = 'Elt isolé???????????????????????????????' WRITE(IOIMP,*) MOTERR(1:40) INDER = 1 GOTO 9999 ENDIF C C*** On parcoure la liste chaîné dans un sens et on la rallonge C de nouveaux voisins dans l'autre C 2 CONTINUE ICEN = KRICEN.LECT(NOCEN1) IF(ICEN .NE. 0)THEN C C******* NOCEN1 est un point centre C IPOS1 = LISTPO.INDEX1(ICEN) NVOIS1 = LISTPO.INDEX1(ICEN+1) - IPOS1 DO IVOIS = 1, NVOIS1 - 1, 1 NOCEN2 = LISTPO.LISPO1(IPOS1+IVOIS) IF(MLEVOI.LECT(NOCEN2) .EQ. 0)THEN C C************** On peut le rajouter dedans MLEVOI.LECT C MLEVOI.LECT(NOCEN2)=LAST0 LAST0 = NOCEN2 ICOORX = ((IDIM + 1) * (NOCEN2 - 1))+1 DIST = (MCOORD.XCOOR(ICOORX) - XCEN)**2 DIST = DIST + ((MCOORD.XCOOR(ICOORX+1)- YCEN)**2) IF(IDIM .EQ. 3) & DIST = DIST + ((MCOORD.XCOOR(ICOORX+2)- ZCEN)**2) C C************* Je cherche la position ou il faut le mettre C DO IVOIS2 = 1, INEUVO, 1 DIST1 = MLRDIS.PROG(IVOIS2) IF(DIST .LT. DIST1) GOTO 3 ENDDO IVOIS2 = INEUVO + 1 3 CONTINUE DO IVOIS3 = INEUVO, IVOIS2, -1 MLEVO2.LECT(IVOIS3+1) = MLEVO2.LECT(IVOIS3) MLRDIS.PROG(IVOIS3+1) = MLRDIS.PROG(IVOIS3) ENDDO MLEVO2.LECT(IVOIS2) = NOCEN2 MLRDIS.PROG(IVOIS2) = DIST INEUVO = INEUVO + 1 ENDIF ENDDO ENDIF NOCEN1 = MLEVOI.LECT(NOCEN1) C C**** Le dernier de la liste chaîné est NOCENE C IF(NOCEN1 .NE. NOCENE) GOTO 2 C C C**** Il faut nettoyer MLEVOI C NOCEN1 = LAST0 5 CONTINUE NOCEN2 = MLEVOI.LECT(NOCEN1) MLEVOI.LECT(NOCEN1) = 0 NOCEN1 = NOCEN2 IF(NOCEN1 .NE. -1) GOTO 5 C IF(INEUVO .LT. (NVOIMI - NTOTP))THEN MOTERR(1:40) = 'Pas assez de noueds voisins ' WRITE(IOIMP,*) MOTERR(1:40) INDER = 1 GOTO 9999 ENDIF C C**** Test 1: la liste ordonnée C C write(*,*) ' ' C write(*,*) 'NOCENE =', NOCENE C write(*,*) 'INEUVO =', INEUVO C do ivois3 = 1, ineuvo, 1 C write(*,*) mlevo2.lect(ivois3), mlrdis.prog(ivois3) C enddo C C**** Maintenant il faut selectioner les noeuds au ajouter C C IPOS1 = MAX(NVOIMI - NTOTP +1,1) IF(IPOS1 .GT. INEUVO)THEN MOTERR(1:40) = 'Pas assez de noueds voisins ' WRITE(IOIMP,*) MOTERR(1:40) INDER = 1 GOTO 9999 ENDIF C C**** REPEAT UNTIL IVOIS < (INEUVO+1) .AND. INDER = 1 C IPOS0 = 1 10 CONTINUE DIST = MLRDIS.PROG(IPOS1) * 1.05D0 DO IVOIS = IPOS1+1, INEUVO, 1 DIST1 = MLRDIS.PROG(IVOIS) IF(DIST1 .GT. DIST) GOTO 4 ENDDO IVOIS = INEUVO + 1 C C IVOIS-1 est le dernier qui je prends C 4 CONTINUE C C IPOS0 est le premier qui n'appartient pas C MLEVTO.LECT C DO IVOIS2 = IPOS0, IVOIS - 1, 1 NTOTP = NTOTP + 1 MLEVTO.LECT(NTOTP)=MLEVO2.LECT(IVOIS2) ENDDO IF(IVOIS .EQ. (INEUVO+1))THEN C C MLEVTO.LECT contient déjà tous les noueds C dedans MLEVO2.LECT C CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, & MATCEL,MLRSIG,MLRCEL,COEF,INDER) C IF(INDER .EQ. 1)THEN C C**** Test 2: Tous les voisins C MOTERR(1:40) = 'Probleme dans la inversion de matrices ' WRITE(IOIMP,*) MOTERR(1:40) ENDIF ELSEIF(IVOIS .LT. (INEUVO+1))THEN CALL GRQUA3(CONTOL,NBMAX,NTOTP,NVOIMI,MLEVTO,MATA,MATU,MATV, & MATCEL,MLRSIG,MLRCEL,COEF,INDER) IF(INDER .EQ. 1)THEN IPOS0 = IVOIS IPOS1 = IVOIS GOTO 10 ELSEIF(INDER .EQ. 2)THEN GOTO 9999 ENDIF ELSEIF(IVOIS .GT. (INEUVO+1))THEN C C Erreur dans la subroutine C CALL ERREUR(5) ENDIF C 9999 RETURN END