Télécharger graco2.eso

Retour à la liste

Numérotation des lignes :

  1. C GRACO2 SOURCE CB215821 17/07/25 12:27:31 9516
  2. SUBROUTINE GRACO2(MMATRX)
  3. IMPLICIT INTEGER(I-N)
  4. IMPLICIT REAL*8(A-H,O-Z)
  5. C TANT QUE OOOVAL(1,4) NE MARCHE PAS SUR CRAY
  6. C
  7. C COPIE DE CHOLE POUR FAIRE UN CHOLEVSKI INCOMPLET DE PLUS ON GARDE LA
  8. C MATRICE ASSEMBLEE ( SEMGENT LLIGN)
  9. C
  10. PARAMETER (LPCRAY=10000)
  11. INTEGER OOOVAL,OOOLEN
  12. dimension ittime(4)
  13. C
  14. C **** MISE SOUS FORME A=Lt D L DE LA MATRICE MMATRX
  15. C SEULES LES VALEURS INITIALEMENT NON NULLES SONT CALCULEES
  16. C
  17. -INC CCOPTIO
  18. -INC SMMATRI
  19. -INC CCASSIS
  20. -INC CCHOLE
  21. SEGMENT KIVPO(IIMAX)
  22. SEGMENT KIVLO(IIMAX)
  23. SEGMENT ITMASQ(NBLIG)
  24. SEGMENT IMASQ(LMASQ)
  25.  
  26. external graco5i
  27. SAVE IPASV
  28. DATA IPASV/0/
  29. C character*8 zen
  30. C equivalence (zen,izen)
  31. call timespv(ittime)
  32. kcour=(ittime(1)+ittime(2))/10
  33. kcourp=kcour
  34. kcouri=kcour
  35. kdiff=0
  36. kcour=0
  37. perf=0
  38. perfp=-1
  39. nbchan=1
  40. nbopit=0
  41. iposm=0
  42. C zen='CPU'//char(0)
  43. C le=4
  44. lolig=0
  45. nvaor=0
  46. nbthro=nbthrs
  47. ithrd=0
  48. if (nbthro.gt.1) then
  49. ithrd=1
  50. call threadii
  51. endif
  52. nbthr=nbthro
  53. do ith=1,nbthro
  54. nbop(ith)=0
  55. enddo
  56. C* en fait la parallelisation n'apporte rien car la matrice est trop creuse
  57. C* nbthr=1
  58. MACTIT=OOOVAL(1,1)
  59. prec=1D-30
  60. nvstrm=0
  61. MMATRI=MMATRX
  62. SEGACT,MMATRI*MOD
  63. MILIG1=IILIGN
  64. IASLIG= IILIGN
  65. SEGINI,MILIGN=MILIG1
  66. IILIGN = MILIGN
  67. INO=ILIGN(/1)
  68. MDIA1=IDIAG
  69. IASDIA=IDIAG
  70. SEGINI,MDIAG=MDIA1
  71. C------------------------
  72. IDIAG = MDIAG
  73. C-------------------------
  74. NBLIG=INO
  75. SEGINI ITMASQ
  76. C* comme la matrice est beaucoup plus creuse qu'en direct, on prend un masdim plus petit
  77. precc=prec
  78. INC=DIAG(/1)
  79. INCC=INC
  80. MIMIK=IIMIK
  81. MINCPO=IINCPO
  82. SEGACT,MINCPO,MIMIK
  83. IPLUMI=IMIK(/2)*2 +4
  84. IL2=0
  85. LTRK=MAX(LPCRAY+1,OOOVAL(1,4))
  86. IIMAX=((((IJMAX+IPLUMI) +LTRK)/LTRK)+1)*LTRK
  87. SEGINI KIVPO,KIVLO
  88. INEG=0
  89. NBLAG=0
  90. NENSLX=0
  91. NVSTOC=0
  92. NVSTOR=0
  93. C
  94. C ngmaxy vient de option
  95. C ngmaxx est la valeur autoajustee
  96. C ngmaxz est la valeur utilisee dans les tests (grande si il y a du debordement)
  97. ngmaxx=ngmaxy
  98. NGMAXZ=NGMAXX
  99. C
  100. C
  101. C **** DEBUT DE LA TRIANGULARISATION. ON PREND NOEUD A NOEUD,
  102. C **** DECOMPACTAGE PUIS TRAVAIL SUR LES LIGNES DU NOEUDS
  103. C
  104. C **** LA LONGUEUR DE LA PLUS GRANDE LIGNE EST DONNEE PAR IMAX
  105. C
  106. 1 CONTINUE
  107. IVALMA=IJMAX+IPLUMI
  108. IL1=IL2+1
  109. IMIN=IL1
  110.  
  111.  
  112. mactic=0
  113. DO 2 I=IL1,INO
  114.  
  115. LLIGN= ILIGN(I)
  116. SEGACT /ERR=3/LLIGN
  117.  
  118.  
  119. NA= IMMMM(/1)
  120. C* write (6,*) ' chole ligne noeud inconnues ',i,ipno(i),na
  121. NBPAR=NA+1
  122. NVALL= NJMAX
  123. lolig=nvall/na+1
  124. if (lolig.gt.6000.and.lolig.gt.2*loligp.and.
  125. > mactic.gt.NGMAXZ*nbthro/2 .and.
  126. > i-il1.ge.nbthro) goto 3
  127. loligp=lolig
  128.  
  129. SEGINI /ERR=3/ LIGN
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. C recuperer la longueur du segment
  137. mactic=mactic+ooolen(lign)
  138. C* if ((mactic.gt.mactit/2.or.mactic.gt.mactit-nvstrm*12)
  139. C* > .AND.I.GT.IL1) THEN
  140. C* SEGSUP LIGN
  141. C* GOTO 3
  142. C* endif
  143. NVSTOC=NVSTOC + NVALL
  144. IVALMA=IVALMA + NVALL
  145.  
  146.  
  147. nvaor = nvaor + xxva(/1)
  148. C
  149. C **** DECOMPACTAGE
  150. C
  151. LMASQ=NVALL/MASDIM+2
  152. SEGINI IMASQ
  153. IPA=1
  154. valmi=1.D50
  155. valma=1.D-50
  156. DO 121 JPA=1,NA
  157. IVPO(JPA)=IPA
  158. KPA=IPPO(JPA+1)-IPPO(JPA)
  159. IPP=IPPO(JPA)
  160. IPPVV(JPA)=IPA-1
  161. LPA=LDEB(JPA)
  162. LPA1=LPA-IPA
  163. DO 122 MPA=1,KPA
  164. LL=LINC(MPA+IPP)
  165. IPLA=LL-LPA1
  166. VAL(IPLA)=XXVA(MPA+IPP)
  167. IMASQ(IPLA/MASDIM+1)=1
  168.  
  169. XVAL1 = abs(val(ipla))
  170. valma=max(valma,XVAL1)
  171. if (val(ipla).ne.0.D0) valmi=min(valmi,XVAL1)
  172. 122 CONTINUE
  173. C REDUCTION EVENTUELLE DE LA LIGNE A SA PARTIE SIGNIFICATIVE
  174. C > 1E-50*DIAGONALE
  175. C*** ***
  176. VALREF=ABS(VAL(IPA+IMMMM(JPA)-LPA))*1D-50
  177. ITROP=0
  178. IPAHAU=IPA+IMMMM(JPA)-LPA
  179. C* DO 123 ILD=IPA,IPAHAU-1
  180. C* IF (ABS(VAL(ILD)).GT.VALREF) GOTO 124
  181. C* ITROP=ITROP+1
  182. C*23 CONTINUE
  183. C124 CONTINUE
  184. C* ITROP=0
  185. C* IF (ITROP.NE.0) THEN
  186. C* LPA=LPA+ITROP
  187. C* JILD=IPAHAU-ITROP
  188. C* DO 125 ILD=IPA,JILD
  189. C* VAL(ILD)=VAL(ILD+ITROP)
  190. C125 CONTINUE
  191. C* DO 126 ILD=JILD+1,JILD+ITROP
  192. C* VAL(ILD)=0.D0
  193. C126 CONTINUE
  194. C* NVALL=NVALL-ITROP
  195. C* NVSTOC=NVSTOC-ITROP
  196. C* IVALMA=IVALMA-ITROP
  197. C* ENDIF
  198. IPA=IPA+IMMMM(JPA)-LPA+1
  199. IMMM(JPA)=IPNO(LPA)
  200. IF(IMIN.GT.IPNO(LPA)) IMIN=IPNO(LPA)
  201. 121 CONTINUE
  202. C if (valmi.le.1d-20) then
  203. C write (6,*) ' valmi valma diag ',valmi,valma,val(val(/1))
  204. C endif
  205. IF (NVALL.NE.VAL(/1))SEGADJ LIGN
  206. C*** **** ****
  207. C* LMASQ=NVALL/MASDIM+1
  208. C* SEGINI IMASQ
  209. C* DO ILD=1,VAL(/1)
  210. C* IF (ABS(VAL(ILD)).GT.VALREF) IMASQ(ILD/MASDIM+1)=1
  211. C* ENDDO
  212. ITMASQ(I)=IMASQ
  213.  
  214. if (na.gt.0) then
  215. IPREL= IMMMM(1)
  216. IDERL= IMMMM(NA)
  217. endif
  218. IPPVV(NA+1)=IPA-1
  219. C SEGSUP LLIGN
  220. C SEGDES LLIGN
  221. ILIGN(I)=LIGN
  222. IF(IIMPI.EQ.1525) THEN
  223. WRITE( IOIMP,4987) I
  224. 4987 FORMAT (' NOEUD NUMERO ',I5)
  225. LL=VAL(/1)
  226. WRITE(IOIMP, 4926)(VAL(MPA),MPA=1,LL)
  227. 4926 FORMAT(' VAL ' , 10E11.4)
  228. ENDIF
  229. nbths=nbthr
  230. if (mod(i+1-il1,nbths).eq.0
  231. > .and.mactic.gt.NGMAXZ*nbthro.OR.I+1-il1.gt.1000) then
  232. il2=i
  233. GOTO 4
  234. endif
  235. 2 CONTINUE
  236. IL2=INO
  237. GO TO 4
  238. 3 IL2=I-1
  239. 4 CONTINUE
  240. nbthro=nbthrs
  241. if (mactic.le.50000) nbthro=1
  242. nbthr=nbthro
  243.  
  244. IF(IL2.GE.IL1) GO TO 40
  245. C
  246. C **** APPEL AUX ERREURS MESSAGE PAS ASSEZ DE PLACE MEMOIRE
  247. C
  248. ITYP=48
  249. CALL ERREUR(48)
  250. if (ithrd.eq.1) call threadis
  251. RETURN
  252. 40 CONTINUE
  253. IM=INC
  254.  
  255.  
  256.  
  257.  
  258. DO 352 IH=IL2 ,IL1,-1
  259. LIGN= ILIGN(IH )
  260. IL=INC
  261. DO 354 JH=1, IMMM(/1)
  262. IM=MIN(IM, IMMM(JH))
  263. IL=MIN(IL, IMMM(JH))
  264. 354 CONTINUE
  265. IML=IL
  266. IMM=ipno(IM)
  267. 352 CONTINUE
  268. C 353 CONTINUE
  269. LIGN=ILIGN(IL1)
  270. IL11=IPREL
  271. C
  272. C **** BOUCLE *5* TRAVAILLE SUR LE NOEUD I QUI EST EN LECTURE
  273. C
  274. C if (ngmaxz.eq.ngmaxy) write (6,*) 'resolution in core ',il1,il2
  275. C if (ngmaxz.ne.ngmaxy) write (6,*) 'resolution out of core ',
  276. C > il1,il2
  277. ipos=0
  278. iper=imin
  279. ider=imin-1
  280. ngmaxz=ngmaxx
  281. macsec=0
  282. isec=0
  283.  
  284. DO 5 I=IMIN,IL2
  285. if (ierr.ne.0) return
  286. LIG1=ILIGN(I)
  287. IMASQ=ITMASQ(I)
  288. IF(I.LT.IL1) GO TO 7
  289. C
  290. C ******* LE NOEUD I EST EN MEMOIRE IL EST TRIANGULE JUSQU'A
  291. C ******* IPREL IL FAUT CONTINUER TOUTE LES LIGNES PUIS CALCULER
  292. C ******* LE TERME DIAGONAL
  293. C
  294. LIGN=LIG1
  295. DO 156 KHG=1,IMMM(/1)
  296. II=IPREL-1+KHG
  297. IMMM(KHG)=0
  298. NN=IPPVV(KHG+1)
  299. NNM1=IPPVV(KHG)
  300. N=NN-NNM1
  301. DIAG(II)=VAL(NN)
  302. IF(N.EQ.1) GO TO 8
  303. NMI=N-II
  304. IDEP=MAX(IL11,2-NMI)
  305. KIDEP=IDEP+NMI
  306. KI1=N-1
  307. KQ=-NMI
  308. VAL(NN)=VAL(NN)+
  309. # GRACO3(ILIGN,LIGN,VAL(1+IPPVV(KHG)),DIAG(1-NMI),IPNO(1-NMI),
  310. # IPPVV(1),KHG,IVPO(1),KIDEP,KI1,KQ,IMASQ(1),1+IPPVV(KHG),
  311. # PREC,ittr(1),nc)
  312. imasq(nn/masdim+1)=1
  313. imasq(n/masdim+1)=1
  314. 8 CONTINUE
  315. if (diag(ii).ne.0.d0) then
  316. if (val(nn).lt.abs(diag(ii))*0.5) val(nn)=diag(ii)
  317. endif
  318.  
  319. IF(ITTR(II).EQ.0.AND.
  320. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-10) GO TO 12
  321. IF(ITTR(II).NE.0.AND.
  322. & ABS(VAL(NN)).GT.ABS(DIAG(II))*1.E-6) GO TO 12
  323. C
  324. C **** ON VIENT DE DETECTER UN MODE D'ENSEMBLE
  325. C **** ON AJOUTE A LA STRUCTURE UN RESSORT EGAL A CELUI QUI EXISTAIT
  326. C **** AU PREALABLE SUR CETTE INCONNUE.
  327. C
  328. VAL(NN)=DIAG(II)*1d-4
  329. NENS=NENS+1
  330. IMMM(KHG)=NENS
  331. if(ittr(ii).NE.0) NENSLX=NENSLX+1
  332. 12 CONTINUE
  333. IMASQ(NN/MASDIM+1)=1
  334. DIAG(II)=VAL(NN)
  335. IF(DIAG(II).NE.0.d0) GO TO 41
  336. KQ1=1+NNM1
  337. KQN=N+NNM1
  338. DO 16 LFG=KQ1,KQN
  339. if (imasq(lfg/masdim+1).le.0) goto 16
  340. IF(VAL(LFG).NE.0.d0) GO TO 17
  341. 16 CONTINUE
  342. DIAG(II)=1.D0
  343. if (ittr(ii).ne.0) diag(ii)=-1d0
  344. VAL(NN)=DIAG(II)
  345. GO TO 41
  346. 17 CONTINUE
  347. C
  348. C **** ENVOI ERREUR MATRICE SINGUIERE
  349. C
  350. ITYP=49
  351. INTERR(1)=I
  352. CALL ERREUR(49)
  353. if (ithrd.eq.1) call threadis
  354. RETURN
  355. C
  356. C **** ON COMPTE LE NOMBRE DE TERMES DIAGONAUX NEGATIFS
  357. C ET LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  358. C
  359. 41 IF(DIAG(II).LT.0.D0) INEG=INEG+1
  360. IF(ITTR(II).NE.0) NBLAG=NBLAG+1
  361. 156 CONTINUE
  362. NA=IMMM(/1)
  363. C
  364. C RECOMPACTAGE DE LIGN (DEJA ENTIEREMENT TRAITEE)
  365. C
  366. izro=2
  367. CALL COMPAC(VAL(1),NBPAR,KIVPO(1),KIVLO(1),
  368. # NVALL,IPPVV(1),izro ,NA ,1D-30,imasq(1),iprel,iderl)
  369. C on recree lign car la compacter en place emiette la memoire
  370. segadj lign
  371. segini,lig1=lign
  372. segsup lign
  373. lign=lig1
  374. ilign(i)=lign
  375.  
  376. NVSTOR=NVSTOR+NVALL
  377. nvstrm=max(nvstrm,nvall)
  378. DO 143 LHG=1,NBPAR
  379. IVPO(2*LHG-1)=KIVPO(LHG)
  380. IVPO(2*LHG)=KIVLO(LHG)
  381. 143 CONTINUE
  382. segsup imasq
  383. imasq=itmasq(i)
  384. if (i.gt.1) then
  385. lig1=ilign(i-1)
  386. segdes lig1
  387. iderac=min(iderac,i-2)
  388. endif
  389. C* WRITE (6,*) ' NOEU ',I,' DIAG INVERSE ',VAL(VAL(/1))
  390. C* WRITE (6,*) ' LIGNE APRES COMPACTION VAL '
  391. C* WRITE (6,*) (VAL(IK),IK=1,VAL(/1))
  392. C* WRITE (6,*) 'KIVPO '
  393. C* WRITE (6,*) (KIVPO(IK),IK=1,NBPAR)
  394. C* WRITE (6,*) 'KIVLO '
  395. C* WRITE (6,*) (KIVLO(IK),IK=1,NBPAR)
  396.  
  397. C **** ON TRIANGULARISE LES AUTRES LIGNES
  398. C
  399. il1=il1+1
  400. if (il1.gt.il2) goto 5
  401. LIG1=ILIGN(I)
  402. lign=ilign(il1)
  403. IL11=IPREL
  404. goto 7
  405. C 72 continue
  406. 71 continue
  407. if (i.ne.il1-1) then
  408. ngmaxz= (oooval(1,1)/nbths)*0.66667
  409. C** write (6,*) ' passage au grand ngmax ',i,il1,il2,ngmaxz
  410. endif
  411. if (isec.ne.0) ider=isec
  412. macsec=0
  413. isec=0
  414. C soit parce qu'on a fini, soit parce qu'on manque de memoire
  415. C il faut executer les lignes activees puis les desactiver
  416. C lancer les chole3 et attendre qu'ils soient finis
  417. if (ipos.ne.0) then
  418. C write (6,*) ' lancement thread ',iper,ider,il1,il2
  419. if (iper.gt.ider) then
  420. write (6,*) ' erreur interne chole '
  421. write (6,*) ' noeud ',i,' sur ',ino
  422. write (6,*) ' nombre de termes actuels ',nvstor
  423. call erreur(5)
  424. endif
  425. if (nbthr.gt.1) then
  426. do ith=2,nbthr
  427. call threadid(ith,graco5i)
  428. enddo
  429. call graco5i(1)
  430. do ith=2,nbthr
  431. call threadif(ith)
  432. enddo
  433. C en multithread il peut y avoir n'importe quoi dans oov(1) du
  434. C aux acces simultanes et ca crache gemat. donc :
  435. oov(1)=0
  436. else
  437. do ith=1,nbthr
  438. call graco5i(ith)
  439. enddo
  440. endif
  441. endif
  442.  
  443. C test ctrlC
  444. if (ierr.ne.0) goto 9999
  445. iposm=max(iposm,ipos)
  446. ipos=0
  447. iderf=ider-1
  448. if (ider.ne.il1-1) iderf=ider
  449. do il=iderf,iper,-1
  450. lign=ilign(il)
  451. segdes lign
  452. C write (6,*) ' desactivation ligne ',il
  453. enddo
  454. iderac=min(iderac,iper-1)
  455. iper=ider+1
  456. C write (6,*) ' iper ider il1 ',iper,ider,il1
  457. if (iper.ne.il1) goto 7
  458.  
  459. goto 5
  460. 7 CONTINUE
  461. C blocage secondaire sur les lignes en lecture pour minimiser le cache
  462. if (isec.le.iper+nbthrs) then
  463. if (macsec.gt.ngmaxz*nbthrs*2) then
  464. isec=i-1
  465. C on continue a activer pour voir si on doit passer en mode lent
  466. ** goto 72
  467. endif
  468. endif
  469. C* if (i.gt.iderac) SEGACT/err=71/LIG1
  470. SEGACT/err=71/LIG1
  471. C macsec=macsec+ooolen(lig1)
  472. ipos=ipos+1
  473. ider=i
  474. if (i.gt.iderac) iderac=i
  475. if (i.eq.il1-1) goto 71
  476. 5 CONTINUE
  477. C write (6,*) ' il1 il2 apres 5 ',il1,il2
  478. C DO 11 I=IL1,IL2
  479. C LIGN= ILIGN(I)
  480. SEGDES,LIGN
  481.  
  482.  
  483. C 11 CONTINUE
  484. nbopt=0
  485. do ith=1,nbthro
  486. nbopt=nbopt+nbop(ith)
  487. nbop(ith)=0
  488. enddo
  489. nbopin=nbopt
  490. nbopit=nbopit+nbopin
  491. call timespv(ittime)
  492. kcourp=kcour
  493. kcour=(ittime(1)+ittime(2))/10
  494. kdiff=kcour-kcourp
  495. C* write (6,*) ' nb operation temps ',nbopin,kdiff
  496. if (kdiff.gt.5) then
  497. perf=nbopin/kdiff
  498. C* write (6,*) 'perf ngmaxy il1 il2',perf,ngmaxy,il1i,il2
  499. if (perf.lt.perfp*0.90 .and.nbchan.ne.1 ) then
  500. nbchan=1
  501. ngmaxx=ngmaxx*0.90
  502. perfp=perf
  503. elseif (nbchan.eq.0) then
  504. nbchan=-1
  505. ngmaxx=ngmaxx*1.10
  506. perfp=max(perf,perfp)
  507. else
  508. nbchan=0
  509. endif
  510. C* nbchan=0
  511. ngmaxx=max(10000,min(100000000,ngmaxx))
  512. ngmaxz=ngmaxx
  513. endif
  514.  
  515. iderac=min(iderac,il1-1)
  516. macsec=0
  517. isec=0
  518. IF(IL2.LT.INO) GO TO 1
  519. C ON MET A JOUR LE NOMBRE DE TERMES DIAGONAUX NEGATIF
  520. C ON ENLEVE LE NOMBRE DE MULTIPLICATEUR DE LAGRANGE
  521. C INEG=INEG-NBLAG
  522. C on ne compte pas 2 fois les multiplicateurs qui vont etre
  523. C elimines lors de la resolution car mode d'ensemble
  524. INEG=INEG-(NBLAG-NENSLX)
  525. if (iimpi.ne.0.and.NENSLX.gt.0) WRITE(IOIMP,4820) NENSLX
  526. 4820 FORMAT(I12,' MODES D ENSEMBLE PORTANT SUR DES MULTIPLICATEURS',
  527. 1' DE LAGRANGE DETECTES')
  528.  
  529. IF(IIMPI.EQ.1) WRITE(IOIMP,4821) NVSTOC
  530. 4821 FORMAT( ' NOMBRE DE VALEURS DANS LE PROFIL',I12)
  531. IF(IIMPI.EQ.1) WRITE(IOIMP,4822) NVSTOR
  532. 4822 FORMAT( ' NOMBRE DE VALEURS STOCKEES DANS LE PROFIL',I12)
  533. INTERR(1)=NVSTOR
  534. reaerr(1)=nvstor/inc**(4./3)
  535. reaerr(2)=2*nbopit/1D3/max(1,(kcour-kcouri))
  536. IF (IPASV.EQ.0) CALL ERREUR(-278)
  537. IPASV=1
  538. SEGDES,MINCPO
  539. SEGDES,MIMIK
  540. SEGDES,MMATRI
  541. SEGDES,MILIGN
  542. SEGDES,MDIAG
  543. MMATRX=MMATRI
  544. SEGSUP KIVPO,KIVLO,ITMASQ
  545. 9999 continue
  546. if (ithrd.eq.1) call threadis
  547. RETURN
  548. END
  549.  
  550.  
  551.  
  552.  
  553.  
  554.  
  555.  
  556.  
  557.  
  558.  
  559.  
  560.  
  561.  
  562.  
  563.  
  564.  
  565.  
  566.  
  567.  
  568.  
  569.  
  570.  
  571.  
  572.  
  573.  
  574.  
  575.  
  576.  
  577.  
  578.  
  579.  
  580.  
  581.  
  582.  
  583.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales