C SEPA2 SOURCE PV090527 25/11/03 21:15:15 12394 C**************************************************************************** C**************************************************************************** C***************SEPAR...trouve separation************************************ C**************************************************************************** C**************************************************************************** C SEPAR trouve la separation a partir du domaine defini par C MASQ=.TRUE. et du noeud appele PIVOT, renvoie DIMSEP,le nombre de C noeuds contenant dans la separation, MASQ=.FALSE. pour les noeuds C appartenant a la separation, renumerote celle-ci dans IADJ, SUBROUTINE SEPA2(IADJ,JADJC,PIVOTA,MASQUE,DIMSEP,N,IPOS, & NODES,IPOSMAX,nrelong,noelon,noel2, > londim,nbthr,icco) IMPLICIT INTEGER(I-N) -INC PPARAM -INC CCOPTIO -INC CCREEL common/cnumop/fcout2c(128),fcoutc(128),iadjc,jadjcc,nrelongc, > noelonc,isensc(128),dimlonc(128),masqc,nodesc,iposc,nbtotc, > lfrontc(128),londimc,lpivc(128),ipointc,noel2c, > boolc(128),ldim2c(128),nbthrc,pivotc(50,128),npointc(128), > npoint2c(128),imaxc,lfronc(128),iccon,ncouch integer iths(128) external noepr2i,noepr2j integer dimlonc,boolc SEGMENT IADJ(0) SEGMENT JADJC(0) integer pivot(50),pivot2(50),pivota,pivots,pivotc,pivotb SEGMENT MASQUE LOGICAL MASQ(0) ENDSEGMENT SEGMENT IPOS(0) INTEGER NODES INTEGER BOOL INTEGER DIMSEP,N REAL*8 FCOUT,FCOUT2,fcout2c,fcoutc real*8 xran INTEGER LONG,LONG2,L,DIMLON C LONG2,LONG correspond au pseudo-diametre. C L determine les noeuds appartenant a la separation. SEGMENT NRELONG(NODES) C NRELONG contient pour chaque noeud sa profondeur. SEGMENT NOELON(NODES) SEGMENT NOEL2(NODES) SEGMENT LONDIM(NODES) C NOELON contient les noeuds de profondeur LONG. C DIMLON= dimension de NOELON. C Initialisation des segments de travail. * write (6,*) ' entree dans sepa2 nbthr ',nbthr * noepr2 travaille differamment suivant que nbtot est nul ou pas * et aussi en fonction de iccon FCOUT=0.D0 FCOUT2=0.D0 iccon=icco nbtot=0 lpiv=1 pivot(lpiv)=pivota MDOMN=IPOS(PIVOTA+NODES) CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1), > isens,DIMLON,NODES,IPOS(1+nodes),NBTOT,lfront,lfron, > londim(1),fcout,lpiv,iccon,icouch) do j=1,dimlon noel2(j)=noelon(j) enddo * insertion d'un deuxieme pivot pour faire les partitions par la suite * on part du point extreme et on reprend le nouveau point extreme if (dimlon.ge.2) then lpiv=1 pivot(lpiv)=noelon(dimlon) do ix=1,dimlon ii=noelon(ix) if(ii.ne.0) nrelong(ii)=0 enddo nbtot=0 CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1), > isens,DIMLON,NODES,IPOS(1+nodes),NBTOT,lfront,lfron, > londim(1),fcout,lpiv,iccon,icouch) ncouch=icouch lpiv=2 pivot(lpiv)=noelon(dimlon) endif if (dimlon.eq.0) return * write (6,*) ' sepa2 apres noepr dimlon nbtot fcout ',dimlon, * > nbtot,fcout,pivot(lpiv) * do ix=1,dimlon ii=noelon(ix) if(ii.ne.0) nrelong(ii)=0 enddo if (fcout.lt.2) goto 11 if (dimlon.le.13) goto 11 ldim2=londim(2) * recherche aleatoire * prevoir 2 cases pour les nouveaux pivots car 2 zones npoint=4 do i=lpiv+1,npoint pivot(i)=0 enddo lpiv=npoint bool=2 ipoint=0 imax=0 499 continue do 500 itent=1,igrand if (fcout.lt.2) goto 510 if (bool.le.0 ) goto 510 bool=bool-1 ipoint=ipoint+1 if (ipoint.gt.npoint) then imax=1 ipoint=1 endif npoint2=npoint nbthrl=min(nbthr,nbtot/16+1) ** nbthrl=1 pivots=pivot(ipoint) iadjc=iadj jadjcc=jadjc nrelongc=nrelong noelonc=noelon masqc=masque nodesc=nodes iposc=ipos nbtotc=nbtot londimc=londim ipointc=ipoint noel2c=noel2 nbthrc=nbthrl imaxc=imax do 521 ith=1,nbthrl isensc(ith)=isens dimlonc(ith)=dimlon lfrontc(ith)=lfront lfronc(ith)=lfron fcout2c(ith)=fcout2 lpivc(ith)=lpiv boolc(ith)=bool fcoutc(ith)=fcout ldim2c(ith)=ldim2 npointc(ith)=npoint npoint2c(ith)=npoint2 do jj=1,lfron londim(jj+(ith-1)*nodes)=londim(jj) enddo do jj=1,npoint2 pivotc(jj,ith)=pivot(jj) enddo if (ith.ne.nbthrl.and.dimlon.gt.64 ) then * write (6,*) ' appel threadid ',ith iths(ith)=1 call threadid(ith,noepr2i) else call noepr2i(ith) iths(ith)=0 endif 521 continue isig=2*mod(ipoint,2)-1 do 522 ith=1,nbthrl if (iths(ith).eq.1) call threadif(ith) iths(ith)=0 if (fcoutc(ith).lt.fcout) then if (fcoutc(ith).ne.fcout) bool=2 isens=isensc(ith) fcout=fcoutc(ith) ldim2=ldim2c(ith) npoint=npointc(ith) npoint2=npoint2c(ith) pivot(ipoint)=pivotc(ipoint,ith) pivot(npoint2)=0 endif 522 continue if (ierr.ne.0) return npoint=npoint2 lpiv=npoint 500 continue 510 continue * write (6,*) ' nb pts maitres ',npoint * petit gradient local parallelise + * on garde la supression eventuelle if(dimlon.lt.64) goto 570 isucc=1 ipoins=npoint do 560 itent=1,igrand if (isucc.eq.0) goto 570 if (fcout.lt.2) goto 570 isucc=0 do 550 ipoinu= npoint,1,-1 ipoint=mod(ipoins+ipoinu-2,npoint)+1 lpoint=pivot(ipoint) if (lpoint.le.0) goto 550 iadh=iadj(lpoint+1)-1 iadb=iadj(lpoint) * pour ne pas perdre trop de temps avec un super element iadb=max(iadb,iadh-256) do 580 kb=iadh,iadb,-nbthrl iadjc=iadj jadjcc=jadjc nrelongc=nrelong noelonc=noelon masqc=masque nodesc=nodes iposc=ipos nbtotc=nbtot londimc=londim ipointc=ipoint noel2c=noel2 nbthrc=nbthrl imaxc=imax do 581 ith=nbthrl,1,-1 kk=kb-ith+1 if (kk.lt.iadj(lpoint)) goto 581 if (kk.eq.iadj(lpoint+1)) then pivot(ipoint)=0 else k=jadjc(kk) if (k.eq.0) goto 581 if (IPOS(k+NODES).ne.mdomn) goto 581 pivot(ipoint)=k endif isensc(ith)=isens dimlonc(ith)=dimlon lfrontc(ith)=lfront lfronc(ith)=lfron fcout2c(ith)=fcout2 lpivc(ith)=lpiv boolc(ith)=bool fcoutc(ith)=fcout ldim2c(ith)=ldim2 npointc(ith)=npoint npoint2c(ith)=npoint2 do jj=1,lfron londim(jj+(ith-1)*nodes)=londim(jj) enddo do jj=1,npoint2 pivotc(jj,ith)=pivot(jj) enddo if (ith.ne.1.and.dimlon.gt.64 ) then * write (6,*) ' appel threadid ',ith call threadid(ith,noepr2j) iths(ith)=1 else call noepr2j(ith) iths(ith)=0 endif pivot(ipoint)=lpoint 581 continue 583 continue do 582 ith=nbthrl,1,-1 if (iths(ith).eq.1) then call threadif(ith) iths(ith)=0 endif ** write(6,*) ' apres noepr2j ',fcout,fcoutc(ith) if (fcoutc(ith).lt.fcout) then ** write(6,*) 'amelioration gradient',fcout, ** > fcoutc(ith),pivotc(ipoint,ith) if (fcoutc(ith).ne.fcout) bool=2 isens=isensc(ith) fcout=fcoutc(ith) ldim2=ldim2c(ith) npoint=npointc(ith) npoint2=npoint2c(ith) pivot(ipoint)=pivotc(ipoint,ith) pivot(npoint2)=0 isucc=1 endif 582 continue if (ierr.ne.0) return npoint=npoint2 lpiv=npoint if (isucc.eq.1) then ipoins=ipoint goto 560 endif 580 continue 550 continue 560 continue 570 continue 11 continue nbtot=0 ** write (6,*) ' sepa2 avant noepr2 final ',pivot(1),pivot(2) CALL NOEPR2(IADJ(1),JADJC(1),PIVOT,NRELONG(1),NOELON(1), > isens,DIMLON,NODES,IPOS(1+nodes),NBTOT,lfront,lfron, > londim(1),fcout,lpiv,iccon,icouch) C PIVOT correspond au noeud pseudo-peripherique. C LONG2 correspond au pseudo-diametre. C DIMSEP=0 C MDOMN est le numéro de mon domaine C pour l'instant,aucun noeud n'appartient a la separation. l=lfront C L correspond a la distance moyenne pour aller d'un bout a l'autre C du domaine. C on arrete de separer si LONG < = 5 ou icouch < = 6 C on masque alors tous les noeuds repondant a cette condition. ** write (6,*) ' sepa2 icouch lfront nbtot ',icouch,lfront,nbtot IF(Lfront.le.1.OR.NBTOT.LE.16) THEN C insertion d'une seule zone IPOSV=IPOS(MDOMN+1) C IPOSV est son rang DO I=0,IPOSMAX IF (IPOS(I+1).GT.IPOSV) IPOS(I+1)=IPOS(I+1)+1 ENDDO IPOSMAX=IPOSMAX+1 if (iposmax.gt.nodes-1) then write (6,*) ' 1 iposmax > nodes ',iposmax,nodes ** iposmax=nodes-1 endif C IPOSV est son rang DO 40 IX=1,DIMLON I=NOELON(IX) IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).NE.0) THEN MASQ(I)=.FALSE. IPOS(I+NODES)=IPOSMAX DIMSEP=DIMSEP+1 ENDIF 40 CONTINUE IPOS(IPOSMAX+1-0)=IPOSV+1 * do 41 i=1,nodes * IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).EQ.0) THEN * write (6,*) ' zone non connexe 1 ',mdomn,i * endif *41 continue * CALL CMG2(IADJ,JADJC,PIVOT,LONG,NRELONG,NOELON, * * DIMLON,MASQUE,NODES,IPOS,DIMSEP) GOTO 50 ENDIF C calcul des nouveaux indices de zones IPOSV=IPOS(MDOMN+1) C IPOSV est son rang DO I=0,IPOSMAX IF (IPOS(I+1).GT.IPOSV) IPOS(I+1)=IPOS(I+1)+3 ENDDO IPOSMAX=IPOSMAX+3 if (iposmax.gt.nodes-3) then write (6,*) ' 2 iposmax > modes ' ** iposmax=nodes-3 endif isens=1 IPOS(IPOSMAX+1-2)=IPOSV+1 IPOS(IPOSMAX+1-1)=IPOSV+1+isens IPOS(IPOSMAX+1-0)=IPOSV+4-isens C On calcule la position de chaque noeud I. C si I a une profondeur L, I appartient a la separation. iwr=1 izcti1=igrand izcti2=igrand izcti3=igrand izcti4=igrand izcti5=igrand izcti6=igrand izcti7=igrand izcti8=igrand izcte1=igrand izcte2=igrand izcte3=igrand izcte4=igrand izcte5=igrand izcte6=igrand izcte7=igrand izcte8=igrand * dimseq=dimsep DO 20 IX=1,DIMLON I=NOELON(IX) * L'intérieur IF(MDOMN.EQ.IPOS(I+NODES).AND.(NRELONG(I).GT.0)) THEN IF((NRELONG(I).LT.L)) THEN * recherche de izcti do kk=iadj(i),iadj(i+1)-1 k=jadjc(kk) if (k.ne.0) then if (nrelong(k).eq.0) then ipkn=ipos(ipos(k+nodes)+1) if (izcti1.eq.ipkn) then elseif (izcti1.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=izcti2 izcti2=izcti1 izcti1=ipkn elseif (izcti2.eq.ipkn) then elseif (izcti2.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=izcti2 izcti2=ipkn ELSEIF (IZCTI3.EQ.IPkn) THEN elseif (izcti3.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=ipkn elseif (izcti4.eq.ipkn) then elseif (izcti4.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=ipkn elseif (izcti5.eq.ipkn) then elseif (izcti5.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=ipkn elseif (izcti6.eq.ipkn) then elseif (izcti6.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=ipkn elseif (izcti7.eq.ipkn) then elseif (izcti7.gt.ipkn) then izcti8=izcti7 izcti7=ipkn elseif (izcti8.eq.ipkn) then elseif (izcti8.gt.ipkn) then izcti8=ipkn endif endif endif enddo IPOS(I+NODES)=IPOSMAX-1 IPOS(I+2*NODES)=IPOS(I+2*NODES)+1 * write (6,*) ' zone ',i,IPOS(I+NODES) if (icouch.le.4) masq(i)=.false. GOTO 20 ENDIF IF(NRELONG(I).EQ.L) THEN * pour que le noeud soit bien dans la separation, il faut qu'il soit * connecter a un noeud de profondeur plus élevée (prob des noeuds milieux) DO 100 J=1,IADJ(I+1)-IADJ(I) K=JADJC(IADJ(I)+J-1) C K:voisin de I * if (k.eq.0) write (6,*) ' k peut valoir 0 ' IF(K.EQ.0) GOTO 100 IF (NRELONG(K).EQ.0) GOTO 100 IF(MDOMN.NE.IPOS(k+NODES).and. > ipos(k+nodes).lt.iposmax-2) then write (6,*) ' sepa2 mauvais domaine ',mdomn, > ipos(k+nodes),nrelong(k) GOTO 100 endif IF(NRELONG(K).GT.L) GOTO 110 100 CONTINUE * recherche de izcti do kk=iadj(i),iadj(i+1)-1 k=jadjc(kk) if (k.ne.0) then if (nrelong(k).eq.0) then ipkn=ipos(ipos(k+nodes)+1) if (izcti1.eq.ipkn) then elseif (izcti1.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=izcti2 izcti2=izcti1 izcti1=ipkn elseif (izcti2.eq.ipkn) then elseif (izcti2.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=izcti2 izcti2=ipkn ELSEIF (IZCTI3.EQ.ipkn) THEN elseif (izcti3.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=izcti3 izcti3=ipkn elseif (izcti4.eq.ipkn) then elseif (izcti4.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=izcti4 izcti4=ipkn elseif (izcti5.eq.ipkn) then elseif (izcti5.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=izcti5 izcti5=ipkn elseif (izcti6.eq.ipkn) then elseif (izcti6.gt.ipkn) then izcti8=izcti7 izcti7=izcti6 izcti6=ipkn elseif (izcti7.eq.ipkn) then elseif (izcti7.gt.ipkn) then izcti8=izcti7 izcti7=ipkn elseif (izcti8.eq.ipkn) then elseif (izcti8.gt.ipkn) then izcti8=ipkn endif endif endif enddo IPOS(I+NODES)=IPOSMAX-1 IPOS(I+2*NODES)=IPOS(I+2*NODES)+1 * write (6,*) ' zone ',i,IPOS(I+NODES) if (icouch.le.4) masq(i)=.false. GOTO 20 ENDIF * * L'extérieur IF(NRELONG(I).GT.L) THEN * recherche de izcte do kk=iadj(i),iadj(i+1)-1 k=jadjc(kk) if (k.ne.0) then if (nrelong(k).eq.0) then ipkn=ipos(ipos(k+nodes)+1) if (izcte1.eq.ipkn) then elseif (izcte1.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=izcte5 izcte5=izcte4 izcte4=izcte3 izcte3=izcte2 izcte2=izcte1 izcte1=ipkn elseif (izcte2.eq.ipkn) then elseif (izcte2.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=izcte5 izcte5=izcte4 izcte4=izcte3 izcte3=izcte2 izcte2=ipkn elseif (izcte3.eq.ipkn) then elseif (izcte3.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=izcte5 izcte5=izcte4 izcte4=izcte3 izcte3=ipkn elseif (izcte4.eq.ipkn) then elseif (izcte4.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=izcte5 izcte5=izcte4 izcte4=ipkn elseif (izcte5.eq.ipkn) then elseif (izcte5.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=izcte5 izcte5=ipkn elseif (izcte6.eq.ipkn) then elseif (izcte6.gt.ipkn) then izcte8=izcte7 izcte7=izcte6 izcte6=ipkn elseif (izcte7.eq.ipkn) then elseif (izcte7.gt.ipkn) then izcte8=izcte7 izcte7=ipkn elseif (izcte8.eq.ipkn) then elseif (izcte8.gt.ipkn) then izcte8=ipkn endif endif endif enddo IPOS(I+NODES)=IPOSMAX IPOS(I+2*NODES)=IPOS(I+2*NODES)+1 * write (6,*) ' zone ',i,IPOS(I+NODES) if (icouch.le.4) masq(i)=.false. GOTO 20 ENDIF GOTO 21 * la frontiere 110 CONTINUE * write (6,*) ' zone ',i,IPOS(I+NODES) MASQ(I)=.FALSE. IPOS(I+NODES)=IPOSMAX-2 DIMSEP=DIMSEP+1 GOTO 20 ENDIF * mars 21 CONTINUE 20 CONTINUE * do 23 i=1,nodes * IF(MDOMN.EQ.IPOS(I+NODES).AND.NRELONG(I).EQ.0) THEN * write (6,*) ' zone non connexe 2 ',mdomn,i * endif * 23 continue * write (6,*) ' taille de frontière ',dimsep-dimseq * write (6,*) ' rang des zones ',izcti,izcte ipos1=max(IPOS(IPOSMAX+1-1),IPOS(IPOSMAX+1+0)) ipos2=min(IPOS(IPOSMAX+1-1),IPOS(IPOSMAX+1+0)) if (izcti1.gt.izcte1) then * write(6,*) 'sepa2 ordre 1 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti1.lt.izcte1) then * write(6,*) 'sepa2 ordre 1 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti2.gt.izcte2) then * write(6,*) 'sepa2 ordre 2 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti2.lt.izcte2) then * write(6,*) 'sepa2 ordre 2 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti3.gt.izcte3) then * write(6,*) 'sepa2 ordre 3 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti3.lt.izcte3) then * write(6,*) 'sepa2 ordre 3 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti4.gt.izcte4) then * write(6,*) 'sepa2 ordre 4 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti4.lt.izcte4) then * write(6,*) 'sepa2 ordre 4 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti5.gt.izcte5) then * write(6,*) 'sepa2 ordre 5 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti5.lt.izcte5) then * write(6,*) 'sepa2 ordre 5 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti6.gt.izcte6) then * write(6,*) 'sepa2 ordre 6 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti6.lt.izcte6) then * write(6,*) 'sepa2 ordre 6 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti7.gt.izcte7) then * write(6,*) 'sepa2 ordre 7 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti7.lt.izcte7) then * write(6,*) 'sepa2 ordre 7 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 elseif (izcti8.gt.izcte8) then * write(6,*) 'sepa2 ordre 8 +' IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 elseif (izcti8.le.izcte8) then * write(6,*) 'sepa2 ordre 8 -' IPOS(IPOSMAX+1-1)=ipos2 IPOS(IPOSMAX+1+0)=ipos1 else * write (6,*) ' ordonnancement impossible' * > izcti1,izcti2,izcti3,izcti4,izcti5,izcti6, * > izcte1,izcte2,izcte3,izcte4,izcte5,izcte6 IPOS(IPOSMAX+1-1)=ipos1 IPOS(IPOSMAX+1+0)=ipos2 endif 50 continue do ix=1,dimlon ii=noelon(ix) if(ii.ne.0) nrelong(ii)=0 enddo RETURN END