Télécharger sormat.eso

Retour à la liste

Numérotation des lignes :

  1. C SORMAT SOURCE PV 20/03/30 21:24:33 10567
  2. ************************************************************************
  3. * NOM : sormat.eso
  4. * DESCRIPTION : Sortie d'objets de type RIGIDITE/CHPOINT définissant
  5. * un problème physique sous forme matricielle
  6. * REFERENCES : - The Matrix Market Exchange Formats: Initial Design,
  7. * Boisvert R. F., Pozo R., Remington K. A. (Dec 1996)
  8. * - The Rutherford-Boeing Sparse Matrix Collection,
  9. * Duff I. S., Grimes R. G., Lewis G. L. (Sep 1997)
  10. ************************************************************************
  11. * HISTORIQUE : 7/06/2012 : JCARDO : création de la subroutine
  12. * HISTORIQUE : 4/12/2012 : JCARDO : ajout de la sortie RB
  13. * + formes ELEMEN et TRIANG
  14. * + mots-clés RHS, SOL et RES
  15. * HISTORIQUE : 16/07/2019 : GOUNAND : implem RESU, FORC (MM assemblé)
  16. * HISTORIQUE : 8/11/2019 : JCARDO : COERIG n'etait pas pris en compte
  17. ************************************************************************
  18. * ENTRÉES :: aucune
  19. * SORTIES :: aucune (sur fichier uniquement)
  20. ************************************************************************
  21. * SYNTAXE (GIBIANE) :
  22. *
  23. * SORT 'MAT' MOT1 MOT2 RIG1
  24. * ('TITR' MOT3)
  25. * ('INCO')
  26. * ('GEOM')
  27. * ('FORC' CHP1)
  28. * ('CONN' CHP2)
  29. * ('RESU' CHP3)
  30. * ('SOLU' CHP4) ;
  31. *
  32. ************************************************************************
  33. SUBROUTINE SORMAT
  34.  
  35. IMPLICIT INTEGER(I-N)
  36. IMPLICIT REAL*8(A-H,O-Z)
  37. LOGICAL LDBG
  38. EXTERNAL LONG
  39.  
  40.  
  41. -INC PPARAM
  42. -INC CCOPTIO
  43. *
  44. -INC SMLMOTS
  45. -INC SMCHPOI
  46. -INC SMRIGID
  47. -INC SMMATRI
  48. -INC SMVECTD
  49. POINTEUR ISMBR.MVECTD
  50. POINTEUR INCX.MVECTD
  51. POINTEUR IR.MVECTD
  52. -INC SMELEME
  53. -INC SMCOORD
  54.  
  55. *
  56. * NPOTOT = nombre de noeuds distincts au total
  57. * NELTOT = nombre d'éléments au total
  58. * NINTOT = nombre d'indices au total
  59. * NVATOT = nombre de valeurs au total dans la matrice creuse
  60. * NBCOPR = nombre de composantes primales
  61. * NBCODU = nombre de composantes duales
  62. * NBPTOT = nombre d'inconnues primales au total
  63. * NBDTOT = nombre d'inconnues duales au total
  64. *
  65. *
  66. * DECLARATIONS DES SEGMENTS ET DES POINTEURS
  67. * ------------------------------------------
  68.  
  69. * Stockage morse provenant de l'include SMMATRIK
  70. SEGMENT PMORS
  71. INTEGER IA(NTT+1)
  72. INTEGER JA(NJA)
  73. ENDSEGMENT
  74. SEGMENT IZA
  75. REAL*8 A(NBVA)
  76. ENDSEGMENT
  77. POINTEUR KMORS.PMORS,KMOR2.PMORS
  78. POINTEUR KIZA.IZA,KIZA2.IZA
  79.  
  80. * Tableaux temporaires permettant de concatener des données avant
  81. * de les imprimer en un seul bloc
  82. SEGMENT MIDATA
  83. INTEGER IWRIT(NBWR)
  84. ENDSEGMENT
  85. SEGMENT MRDATA
  86. REAL*8 XWRIT(NBWR)
  87. ENDSEGMENT
  88.  
  89. * Correspondance local/global pour les noeuds
  90. SEGMENT IPOG2L(NPOMAX)
  91.  
  92. * Tableaux temporaires permettant de mémoriser les indices des
  93. * composantes primales/duales de chaque sous-matrice dans la liste
  94. * globale IMIK/IDUA
  95. SEGMENT MCONUM(NRIGEL)
  96. SEGMENT ICONUM
  97. INTEGER IICOPR(NLIGRP),IICODU(NLIGRD)
  98. ENDSEGMENT
  99.  
  100. * Tableaux indiquant le numéro du d.d.l. primal/dual en fonction du
  101. * numéro du noeud et du numéro de la composante
  102. POINTEUR IDDLPR.MINCPO,IDDLDU.MINCPO
  103.  
  104. * Maillage répertoriant tous les noeuds supportant la matrice
  105. POINTEUR MNOEUD.MELEME
  106.  
  107.  
  108. * LISTES DE MOTS-CLÉS
  109. * -------------------
  110. PARAMETER (NCLE=7)
  111. CHARACTER*4 LCLE(NCLE)
  112. DATA LCLE/'TITR','GEOM','INCO','FORC','CONN','SOLU','RESU'/
  113.  
  114. PARAMETER (NTYP=2)
  115. CHARACTER*4 LTYP(NTYP)
  116. DATA LTYP/'ELEM','ASSE'/
  117.  
  118. CHARACTER*20 MYFMT
  119.  
  120. CHARACTER*16 PTRFMT,INDFMT
  121. CHARACTER*20 VALFMT
  122. INTEGER PTRCRD,INDCRD,VALCRD,TOTCRD
  123. PARAMETER (PTRFMT='(6I12)')
  124. PARAMETER (INDFMT='(6I12)')
  125. PARAMETER (VALFMT='(3E25.16)')
  126. PARAMETER (NPTRFMT=5)
  127. PARAMETER (NINDFMT=5)
  128. PARAMETER (NVALFMT=3)
  129.  
  130. * VARIABLES BOOLÉENNES
  131. * --------------------
  132. LOGICAL ZTITR,ZGEOM,ZINCO,ZFORC,ZELEM,ZSOLU,ZRESU,ZNZER
  133. LOGICAL ZOPEN
  134.  
  135.  
  136. * AUTRES DÉCLARATIONS
  137. * -------------------
  138.  
  139. * Format de sortie des fichiers
  140. CHARACTER*17 CEXTN
  141.  
  142. * Chaîne imprimée dans la section TITRE
  143. CHARACTER*128 CTITR,CTIT2
  144. DATA CTITR/' '/
  145.  
  146. * Nom du fichier
  147. CHARACTER*256 NOMFI2
  148.  
  149. * Coefficient multiplicateur d'une sous-rigidite
  150. REAL*8 COEF
  151.  
  152.  
  153. IDIM1=IDIM+1
  154. LDBG=.FALSE.
  155.  
  156. * +---------------------------------------------------------------+
  157. * | |
  158. * | L E C T U R E D E S A R G U M E N T S |
  159. * | |
  160. * +---------------------------------------------------------------+
  161.  
  162. * =========================
  163. * FORMAT STANDARD DE SORTIE
  164. * =========================
  165.  
  166. * CEXTN='mm' / IEXTN=1 : format Matrix Market
  167. * CEXTN='rb' / IEXTN=2 : format Rutherford Boeing
  168. CALL LIRCHA(CEXTN,1,LEXTN)
  169. IF (CEXTN(1:LEXTN).EQ.'MM'.OR.
  170. & CEXTN(1:LEXTN).EQ.'MATRIX_MARKET') THEN
  171. CEXTN='mm'
  172. IEXTN=1
  173. ELSEIF (CEXTN(1:LEXTN).EQ.'RB'.OR.
  174. & CEXTN(1:LEXTN).EQ.'RUTHERFORD_BOEING') THEN
  175. CEXTN='rb'
  176. IEXTN=2
  177. ELSE
  178. WRITE(*,*) 'erreur : mm ou rb ?'
  179. CALL ERREUR(21)
  180. RETURN
  181. ENDIF
  182.  
  183.  
  184. * =========================
  185. * TYPE DE MATRICE EN SORTIE
  186. * =========================
  187.  
  188. * ITYP=1 : matrices élémentaires
  189. * ITYP=2 : matrice assemblée
  190. CALL LIRMOT(LTYP,NTYP,ITYP,1)
  191. IF (IERR.NE.0) RETURN
  192.  
  193.  
  194. * ===========================
  195. * LECTURE DE L'OBJET RIGIDITE
  196. * ===========================
  197.  
  198. CALL LIROBJ('RIGIDITE',MRIGID,1,IRETOU)
  199. IF (IERR.NE.0) RETURN
  200.  
  201.  
  202. * =====================
  203. * LECTURE DES MOTS-CLÉS
  204. * =====================
  205.  
  206. ZTITR=.FALSE.
  207. ZGEOM=.FALSE.
  208. ZINCO=.FALSE.
  209. ZFORC=.FALSE.
  210. ISMBR=0
  211. ZELEM=.FALSE.
  212. ZSOLU=.FALSE.
  213. ZRESU=.FALSE.
  214. INCX=0
  215.  
  216. 1 CALL LIRMOT(LCLE,NCLE,ICLE,0)
  217. IF (ICLE.EQ.0) GOTO 100
  218.  
  219. GOTO (10,11,12,13,14,15,16),ICLE
  220.  
  221. * Mot-clé TITR
  222. 10 CONTINUE
  223. CALL LIRCHA(CTITR,1,LTITR)
  224. IF (IERR.NE.0) RETURN
  225. ZTITR=.TRUE.
  226. GOTO 1
  227.  
  228. * Mot-clé GEOM
  229. 11 CONTINUE
  230. ZGEOM=.TRUE.
  231. GOTO 1
  232.  
  233. * Mot-clé INCO
  234. 12 CONTINUE
  235. ZINCO=.TRUE.
  236. GOTO 1
  237.  
  238. * Mot-clé FORC
  239. 13 CONTINUE
  240. CALL LIROBJ('CHPOINT',IFORC,1,IRETOU)
  241. IF (IERR.NE.0) RETURN
  242. ZFORC=.TRUE.
  243. GOTO 1
  244.  
  245. * Mot-clé ELEM
  246. 14 CONTINUE
  247. ZELEM=.TRUE.
  248. GOTO 1
  249.  
  250. * Mot-clé SOLU
  251. 15 CONTINUE
  252. CALL LIROBJ('CHPOINT',ISOLU,1,IRETOU)
  253. IF (IERR.NE.0) RETURN
  254. ZSOLU=.TRUE.
  255. GOTO 1
  256.  
  257. * Mot-clé RESU
  258. 16 CONTINUE
  259. CALL LIROBJ('CHPOINT',IRESU,1,IRETOU)
  260. IF (IERR.NE.0) RETURN
  261. ZRESU=.TRUE.
  262. GOTO 1
  263.  
  264.  
  265.  
  266.  
  267. c WRITE(*,*) 'IEXTN=',IEXTN
  268. c WRITE(*,*) 'ITYP=',ITYP
  269. c WRITE(*,*) 'ZTITR=',ZTITR
  270. c WRITE(*,*) 'ZGEOM=',ZGEOM
  271. c WRITE(*,*) 'ZINCO=',ZINCO
  272. c WRITE(*,*) 'ZFORC=',ZFORC
  273. c WRITE(*,*) 'ZELEM=',ZELEM
  274. c WRITE(*,*) 'ZSOLU=',ZSOLU
  275. c WRITE(*,*) 'ZRESU=',ZRESU
  276.  
  277.  
  278.  
  279.  
  280.  
  281. * +---------------------------------------------------------------+
  282. * | |
  283. * | M A T R I C E = F I C H I E R . M T X |
  284. * | |
  285. * +---------------------------------------------------------------+
  286. *
  287.  
  288. 100 CONTINUE
  289.  
  290.  
  291. * On s'assure que la matrice n'est pas vide
  292. SEGACT,MRIGID*MOD
  293. NRIGEL=IRIGEL(/2)
  294. IF (NRIGEL.EQ.0) THEN
  295. MOTERR(1:8)='RIGIDITE'
  296. CALL ERREUR(1027)
  297. RETURN
  298. ENDIF
  299.  
  300.  
  301.  
  302. * ======================================
  303. * RÉCUPÉRATION DU NOM DE BASE DU FICHIER
  304. * ======================================
  305.  
  306. INQUIRE(UNIT=IOPER,OPENED=ZOPEN)
  307. IF (.NOT.ZOPEN) THEN
  308. CALL ERREUR(-212)
  309. WRITE(IOIMP,*) '(via OPTI "SORT")'
  310. MOTERR(1:8)='.'//CEXTN(1:7)
  311. CALL ERREUR(705)
  312. RETURN
  313. ENDIF
  314.  
  315. INQUIRE(UNIT=IOPER,NAME=NOMFIC)
  316. CLOSE(UNIT=IOPER,STATUS='DELETE')
  317. CALL LENCHA(NOMFIC,LC)
  318.  
  319.  
  320.  
  321. * ========================================
  322. * SORTIE DE MATRICE SOUS FORME ÉLÉMENTAIRE
  323. * ========================================
  324.  
  325. IF (ITYP.EQ.1) THEN
  326.  
  327.  
  328. * ***********************************
  329. * LISTE DES NOEUDS ET DES COMPOSANTES
  330. * ***********************************
  331.  
  332. * On stocke les noeuds rencontrés dans un MELEME
  333. NPOMAX=nbpts
  334. NBSOUS=0
  335. NBREF=0
  336. NBNN=1
  337. NBELEM=NPOMAX
  338. SEGINI,IPOG2L,MNOEUD
  339. MNOEUD.ITYPEL=1
  340.  
  341. * On stocke les noms des composantes primales/duales dans les
  342. * tableaux IMIK et IDUA (analogues à ceux existant dans SMMATRI
  343. * pour les matrices assemblées). Les tableaux ICONUM (un par
  344. * sous-rigidité) mémorisent pour chaque variable primale/duale
  345. * l'indice de la composante correspondante dans IMIK/IDUA.
  346. SEGINI,MIMIK,MIDUA,MCONUM
  347.  
  348. * (boucle sur les sous-rigidités)
  349. NPOTOT=0
  350. NELTOT=0
  351. NVATOT=0
  352. NINTOT=0
  353. DO IRIG=1,NRIGEL
  354.  
  355. * Construction du MELEME
  356. * ----------------------
  357. IPT1=IRIGEL(1,IRIG)
  358. SEGACT,IPT1
  359. NEL1=IPT1.NUM(/2)
  360. NNO1=IPT1.NUM(/1)
  361. NELTOT=NELTOT+NEL1
  362. DO I=1,NEL1
  363. DO J=1,NNO1
  364. IF (IPOG2L(IPT1.NUM(J,I)).EQ.0) THEN
  365. NPOTOT=NPOTOT+1
  366. IPOG2L(IPT1.NUM(J,I))=NPOTOT
  367. MNOEUD.NUM(1,NPOTOT)=IPT1.NUM(J,I)
  368. ENDIF
  369. ENDDO
  370. ENDDO
  371.  
  372.  
  373. * Construction des listes IMIK/IDUA
  374. * ---------------------------------
  375. DESCR=IRIGEL(3,IRIG)
  376. SEGACT,DESCR
  377. NLIGRP=LISINC(/2)
  378. NLIGRD=LISDUA(/2)
  379. NINTOT=NINTOT+((NLIGRP+NLIGRD)*NEL1)
  380. NVATOT=NVATOT+(NLIGRP*NLIGRD*NEL1)
  381. SEGINI,ICONUM
  382. MCONUM(IRIG)=ICONUM
  383.  
  384. * composantes primales
  385. DO 101 K=1,NLIGRP
  386. DO J=1,IMIK(/2)
  387. IF (IMIK(J).EQ.LISINC(K)) THEN
  388. IICOPR(K)=J
  389. GOTO 101
  390. ENDIF
  391. ENDDO
  392. IMIK(**)=LISINC(K)
  393. IICOPR(K)=IMIK(/2)
  394. 101 CONTINUE
  395.  
  396. * composantes duales
  397. DO 102 K=1,NLIGRD
  398. DO J=1,IDUA(/2)
  399. IF (IDUA(J).EQ.LISDUA(K)) THEN
  400. IICODU(K)=J
  401. GOTO 102
  402. ENDIF
  403. ENDDO
  404. IDUA(**)=LISDUA(K)
  405. IICODU(K)=IDUA(/2)
  406. 102 CONTINUE
  407.  
  408. ENDDO
  409. NBELEM=NPOTOT
  410. SEGADJ,MNOEUD
  411.  
  412.  
  413. * *********************************************************
  414. * LISTE DES INCONNUES (= DDL) PRIMALES/DUALES DE LA MATRICE
  415. * *********************************************************
  416.  
  417. NBPTOT=0
  418. NBDTOT=0
  419.  
  420. NNOE=NPOTOT
  421. MAXI=IMIK(/2)
  422. SEGINI,IDDLPR
  423. MAXI=IDUA(/2)
  424. SEGINI,IDDLDU
  425.  
  426. DO IRIG=1,NRIGEL
  427.  
  428. IPT1=IRIGEL(1,IRIG)
  429. DESCR=IRIGEL(3,IRIG)
  430. ICONUM=MCONUM(IRIG)
  431. NEL1=IPT1.NUM(/2)
  432. NBP1=IICOPR(/1)
  433. NBD1=IICODU(/1)
  434. DO I=1,NEL1
  435.  
  436. * inconnues primales
  437. DO K=1,NBP1
  438. ICOP=IICOPR(K)
  439. INOP=IPOG2L(IPT1.NUM(NOELEP(K),I))
  440. IF (IDDLPR.INCPO(ICOP,INOP).EQ.0) THEN
  441. NBPTOT=NBPTOT+1
  442. IDDLPR.INCPO(ICOP,INOP)=NBPTOT
  443. ENDIF
  444. ENDDO
  445.  
  446. * inconnues duales
  447. DO K=1,NBD1
  448. ICOD=IICODU(K)
  449. INOD=IPOG2L(IPT1.NUM(NOELED(K),I))
  450. IF (IDDLDU.INCPO(ICOD,INOD).EQ.0) THEN
  451. NBDTOT=NBDTOT+1
  452. IDDLDU.INCPO(ICOD,INOD)=NBDTOT
  453. ENDIF
  454. ENDDO
  455.  
  456. ENDDO
  457.  
  458. ENDDO
  459.  
  460.  
  461.  
  462. * ********************************
  463. * => Sortie MATRIX_MARKET élémentaire
  464. * ********************************
  465.  
  466. IF (IEXTN.EQ.1) THEN
  467.  
  468. * Création du fichier .mtx.mm
  469. M=NBPTOT
  470. N=NBDTOT
  471. NNZER=NELTOT
  472. CALL OPENMM(NOMFIC,'mtx',CTITR,
  473. & M,N,NNZER,'matrix elemental real general')
  474.  
  475.  
  476. DO IRIG=1,NRIGEL
  477.  
  478. COEF=COERIG(IRIG)
  479. IPT1=IRIGEL(1,IRIG)
  480. SEGACT,IPT1
  481. DESCR=IRIGEL(3,IRIG)
  482. XMATRI=IRIGEL(4,IRIG)
  483. ICONUM=MCONUM(IRIG)
  484. SEGACT,XMATRI
  485.  
  486. DO KEL=1,IPT1.NUM(/2)
  487.  
  488. NPR1=LISINC(/2)
  489. NDU1=LISDUA(/2)
  490. WRITE(IOPER,FMT='(I12,1X,I12)')
  491. & NPR1,NDU1
  492.  
  493. DO KDU=1,NDU1
  494. ICOD=IICODU(KDU)
  495. INOD=IPOG2L(IPT1.NUM(NOELED(KDU),KEL))
  496. WRITE(IOPER,FMT='(I15)')
  497. & IDDLDU.INCPO(ICOD,INOD)
  498. ENDDO
  499.  
  500. DO KPR=1,NPR1
  501. ICOP=IICOPR(KPR)
  502. INOP=IPOG2L(IPT1.NUM(NOELEP(KPR),KEL))
  503. WRITE(IOPER,FMT='(I15)')
  504. & IDDLPR.INCPO(ICOP,INOP)
  505. ENDDO
  506.  
  507. DO KDU=1,NDU1
  508. DO KPR=1,NPR1
  509. WRITE(IOPER,FMT='(E25.16)')
  510. & COEF*RE(KDU,KPR,KEL)
  511. ENDDO
  512. ENDDO
  513.  
  514. ENDDO
  515.  
  516. SEGDES,IPT1,DESCR,XMATRI
  517. SEGSUP,ICONUM
  518.  
  519. ENDDO
  520.  
  521.  
  522. * ************************************
  523. * => Sortie RUTHERFORD_BOEING élémentaire
  524. * ************************************
  525.  
  526. ELSEIF (IEXTN.EQ.2) THEN
  527.  
  528. * Création du fichier .mtx.rb
  529. PTRCRD=2*NELTOT+1
  530. INDCRD=NBPTOT+NBDTOT
  531. VALCRD=NVATOT
  532. TOTCRD=PTRCRD+INDCRD+VALCRD
  533. MVAR=MAX(NBPTOT,NBDTOT)
  534. NELT=NELTOT
  535. NVARIX=NINTOT
  536. NELTVL=NVATOT
  537. CALL OPENRB(NOMFIC,'mtx',CTITR,
  538. & TOTCRD,PTRCRD,INDCRD,VALCRD,
  539. & MVAR,NELT,NVARIX,NELTVL,
  540. & 'rre','(I15)','(I15)','(E25.16)')
  541. IF (IERR.NE.0) RETURN
  542.  
  543.  
  544. * ECRITURE DES POINTEURS
  545. * ----------------------
  546. IPTR=1
  547. WRITE(IOPER,FMT='(I15)')
  548. & IPTR
  549. DO IRIG=1,NRIGEL
  550.  
  551. IPT1=IRIGEL(1,IRIG)
  552. DESCR=IRIGEL(3,IRIG)
  553.  
  554. DO KEL=1,IPT1.NUM(/2)
  555. IPTR=IPTR+LISINC(/2)
  556. WRITE(IOPER,FMT='(I15)')
  557. & IPTR
  558. IPTR=IPTR+LISDUA(/2)
  559. WRITE(IOPER,FMT='(I15)')
  560. & IPTR
  561. ENDDO
  562.  
  563. ENDDO
  564.  
  565. * ECRITURE DES NUMEROS D'INCONNUES
  566. * --------------------------------
  567. DO IRIG=1,NRIGEL
  568.  
  569. COEF=COERIG(IRIG)
  570. IPT1=IRIGEL(1,IRIG)
  571. DESCR=IRIGEL(3,IRIG)
  572. ICONUM=MCONUM(IRIG)
  573.  
  574. DO KEL=1,IPT1.NUM(/2)
  575.  
  576. DO KDU=1,LISDUA(/2)
  577. ICOD=IICODU(KDU)
  578. INOD=IPOG2L(IPT1.NUM(NOELED(KDU),KEL))
  579. WRITE(IOPER,FMT='(I15)')
  580. & IDDLDU.INCPO(ICOD,INOD)
  581. ENDDO
  582.  
  583. DO KPR=1,LISINC(/2)
  584. ICOP=IICOPR(KPR)
  585. INOP=IPOG2L(IPT1.NUM(NOELEP(KPR),KEL))
  586. WRITE(IOPER,FMT='(I15)')
  587. & IDDLPR.INCPO(ICOP,INOP)
  588. ENDDO
  589.  
  590. ENDDO
  591.  
  592. SEGSUP,ICONUM
  593.  
  594. ENDDO
  595.  
  596. * ECRITURE DES VALEURS
  597. * --------------------
  598. DO IRIG=1,NRIGEL
  599.  
  600. IPT1=IRIGEL(1,IRIG)
  601. SEGACT,IPT1
  602. DESCR=IRIGEL(3,IRIG)
  603. XMATRI=IRIGEL(4,IRIG)
  604. SEGACT,XMATRI
  605.  
  606. DO KEL=1,IPT1.NUM(/2)
  607. DO KDU=1,LISDUA(/2)
  608. DO KPR=1,LISINC(/2)
  609. WRITE(IOPER,FMT='(E25.16)')
  610. & COEF*RE(KDU,KPR,KEL)
  611. ENDDO
  612. ENDDO
  613. ENDDO
  614.  
  615. SEGDES,IPT1,DESCR,XMATRI
  616.  
  617. ENDDO
  618.  
  619. ENDIF
  620.  
  621. SEGSUP,IPOG2L,MCONUM
  622.  
  623.  
  624.  
  625. * ======================================
  626. * SORTIE DE MATRICE SOUS FORME ASSEMBLÉE
  627. * ======================================
  628. *
  629.  
  630. ELSEIF (ITYP.EQ.2) THEN
  631.  
  632. IF (ICHOLE.NE.0) ICHOLE=0
  633.  
  634. * Assemblage des matrices élémentaires
  635. * ************************************
  636. *
  637. * SG 2019/07 : a noter que INORMU ne sert pas car normalement c'est
  638. * la normalisation des multiplicateurs de Lagrange et ceux-ci ne sont
  639. * pas gérés par sormat.eso (voir plus loin).
  640. * De même la normalisation AUTO ne marche pas car on utilise
  641. * l'assembleur non symétrique (asns1.eso via kres9)
  642. * Par contre la normalisation faite à la main semble marcher mais
  643. * pas très utile.
  644. *
  645. IF (NORINC.EQ.0.AND.NORIND.EQ.0) THEN
  646. INORMU=0
  647. ELSE
  648. INORMU=1
  649. ENDIF
  650. *dbg write(ioimp,*) 'NORINC,NORIND,INORMU=',NORINC,NORIND,INORMU
  651. *dbg write(ioimp,*) 'Assemblage...'
  652. CALL KRES9(MRIGID,INORMU)
  653. IF (IERR.NE.0) RETURN
  654. *dbg write(ioimp,*) 'Assemblage fini'
  655.  
  656. * Mise sous forme morse = CSR (Compressed Sparse Row)
  657. * ***************************************************
  658. * /!\ Même après transformation en morse, les inconnues restent
  659. * decrites par les segments IINCPO et IDUAPO du SMMATRI
  660. * pointé par ICHOLE dans le MRIGID d'origine.
  661. * En theorie, la matrice qui sort de l'assemblage par KRES9
  662. * est carrée, structurellement symétrique. On forcera alors
  663. * NBDTOT = NBPTOT mais on distinguera tout de meme les
  664. * inconnues primales et duales via IINCPO et IDUAPO.
  665. *dbg write(ioimp,*) 'Passage au format Morse...'
  666. CALL KRES10(MRIGID,KMORS,KIZA)
  667. IF (IERR.NE.0) RETURN
  668. *dbg write(ioimp,*) 'Passage au format Morse fini'
  669. IF (LDBG) THEN
  670. IMPR=3
  671. CALL INFMAT(KMORS,KIZA,IMPR,IRET)
  672. ENDIF
  673. *
  674. SEGACT,KMORS,KIZA
  675. NBPTOT=KMORS.IA(/1)-1
  676. NBDTOT=NBPTOT
  677. * /!\ KIZA.A est rempli de 0 inutiles en fin
  678. * => à ne pas utiliser pour déterminer NVATOT
  679. NVATOT=KMORS.JA(/1)
  680.  
  681. *dbg CALL ECMORS(KMORS,KIZA,4)
  682.  
  683. C
  684. C - Conversion du second membre en MVECTD
  685. C et initialisation du résultat
  686. C Cette partie est reprise de KRES8.ESO
  687. IF (ZFORC) THEN
  688. SEGACT MRIGID
  689. ICHOLX=ICHOLE
  690. C On vérifie que le second membre doit être dans le dual
  691. NOID=1
  692. CALL CHVNS(ICHOLX,IFORC,ISMBR,NOID)
  693. IF (IERR.NE.0) RETURN
  694. C
  695. C Gestion normalisation et Lagrange (repris de MONDES)
  696. C
  697. SEGACT ISMBR*MOD
  698. MMATRI=ICHOLE
  699. SEGACT MMATRI
  700. * SG 2019/07
  701. * Verif que ITTR=0 (en effet, on ne traite pas les multiplicateurs
  702. * de Lagrange donc on n'en veut pas pour l'instant)
  703. * En particulier, si on veut les traiter, il faut les dualiser avant
  704. * assemblage car c'est aujourd'hui un prérequis de l'assemblage de
  705. * RESO
  706. IF (IILIGN.NE.0) THEN
  707. MILIGN=IILIGN
  708. SEGACT MILIGN
  709. DO II=1,ITTR(/1)
  710. if (ITTR(II).NE.0) goto 666
  711. ENDDO
  712. ENDIF
  713. IF (IILIGS.NE.0) THEN
  714. MILIGN=IILIGS
  715. SEGACT MILIGN
  716. DO II=1,ITTR(/1)
  717. if (ITTR(II).NE.0) goto 666
  718. ENDDO
  719. ENDIF
  720. *
  721. IF(IDNORD.GT.0) THEN
  722. MDNO1=IDNORD
  723. ELSE
  724. MDNO1=IDNORM
  725. ENDIF
  726. SEGACT MDNO1
  727. INC=MDNO1.DNOR(/1)
  728. DO 45 I=1,INC
  729. ISMBR.VECTBB(I)=ISMBR.VECTBB(I)*MDNO1.DNOR(I)
  730. 45 CONTINUE
  731. SEGDES MDNO1
  732. SEGDES MMATRI
  733. ENDIF
  734. C
  735. C - Conversion du résultat (inconnue) en MVECTD
  736. C
  737. IF (ZRESU) THEN
  738. SEGACT MRIGID
  739. ICHOLX=ICHOLE
  740. C Changement de noms d'inconnues pour le résultat (passage primal ->
  741. C dual)
  742. CALL ECROBJ('CHPOINT ',IRESU)
  743. CALL MACHI2(2)
  744. CALL LIROBJ('CHPOINT ',IRESU,1,IRET)
  745. IF (IERR.NE.0) RETURN
  746. C On vérifie que le second membre doit être dans le dual
  747. NOID=1
  748. CALL CHVNS(ICHOLX,IRESU,INCX,NOID)
  749. IF (IERR.NE.0) RETURN
  750. * Peu utile SEGSUP,IRESU
  751. C
  752. C Gestion normalisation et Lagrange (repris de MONDES)
  753. C
  754. SEGACT INCX*MOD
  755. MMATRI=ICHOLE
  756. SEGACT MMATRI
  757. * SG 2019/07
  758. * Verif que ITTR=0 (en effet, on ne traite pas les multiplicateurs
  759. * de Lagrange donc on n'en veut pas pour l'instant)
  760. * En particulier, si on veut les traiter, il faut les dualiser avant
  761. * assemblage car c'est aujourd'hui un prérequis de l'assemblage de
  762. * RESO
  763. IF (IILIGN.NE.0) THEN
  764. MILIGN=IILIGN
  765. SEGACT MILIGN
  766. DO II=1,ITTR(/1)
  767. if (ITTR(II).NE.0) goto 666
  768. ENDDO
  769. ENDIF
  770. IF (IILIGS.NE.0) THEN
  771. MILIGN=IILIGS
  772. SEGACT MILIGN
  773. DO II=1,ITTR(/1)
  774. if (ITTR(II).NE.0) goto 666
  775. ENDDO
  776. ENDIF
  777. *
  778. MDNOR=IDNORM
  779. SEGACT MDNOR
  780. INC=DNOR(/1)
  781. DO 35 I=1,INC
  782. INCX.VECTBB(I)=INCX.VECTBB(I)/DNOR(I)
  783. 35 CONTINUE
  784. SEGDES MDNOR
  785. SEGDES MMATRI
  786. ENDIF
  787. * Un petit calcul de résidu pour la route !
  788. IF (ZFORC.AND.ZRESU.AND.LDBG) THEN
  789. IMVEC=2
  790. ith=oothrd
  791. IF(ITH.NE.0)THEN
  792. IMVEC=0
  793. ENDIF
  794. SEGINI,IR=ISMBR
  795. C r(0)=b-Ax
  796. CALL GMOMV(IMVEC,'N',-1.D0,KMORS,KIZA,INCX,1.D0,IR)
  797. RNRM2 = GNRM2(IR)
  798. WRITE(IOIMP,*) 'SORT MAT : ||R||=',RNRM2
  799. SEGSUP IR
  800. ENDIF
  801.  
  802. * On (re)active les segments pour la suite de la subroutine
  803. * *********************************************************
  804. SEGACT,MRIGID
  805. MMATRI=ICHOLE
  806. SEGACT,MMATRI
  807.  
  808.  
  809. * ******************************
  810. * => Sortie MATRIX_MARKET assemblée
  811. * ******************************
  812.  
  813. IF (IEXTN.EQ.1) THEN
  814.  
  815. * Création du fichier .mtx.mm
  816. M=NBDTOT
  817. N=NBPTOT
  818. NNZER=NVATOT
  819. CALL OPENMM(NOMFIC,'mtx',CTITR,
  820. & M,N,NNZER,'matrix coordinate real general')
  821. &
  822.  
  823. * Ecriture de la matrice
  824. NTVA=0
  825. DO I=1,M
  826. NIVA=KMORS.IA(I+1)-KMORS.IA(I)
  827. DO J=1,NIVA
  828. WRITE(IOPER,FMT='(I12,1X,I12,1X,E25.16)')
  829. & I,KMORS.JA(NTVA+J),KIZA.A(NTVA+J)
  830. ENDDO
  831. NTVA=NTVA+NIVA
  832. ENDDO
  833.  
  834.  
  835. * **********************************
  836. * => Sortie RUTHERFORD_BOEING assemblée
  837. * **********************************
  838.  
  839. ELSEIF (IEXTN.EQ.2) THEN
  840.  
  841. * Création du fichier .mtx.rb
  842. M=NBDTOT
  843. NVEC=NBPTOT
  844. NAUXD=NVATOT
  845. c PTRCRD=CEILING((NBPTOT+1)/NPTRFMT)
  846. c INDCRD=CEILING(NVATOT/NINDFMT)
  847. c VALCRD=CEILING(NVATOT/NVALFMT)
  848. PTRCRD=(NBPTOT+1)/NPTRFMT
  849. INDCRD=NVATOT/NINDFMT
  850. VALCRD=NVATOT/NVALFMT
  851. IF (MOD((NBPTOT+1),NPTRFMT).GT.0) PTRCRD=PTRCRD+1
  852. IF (MOD(NVATOT,NINDFMT).GT.0) INDCRD=INDCRD+1
  853. IF (MOD(NVATOT,NVALFMT).GT.0) VALCRD=VALCRD+1
  854. TOTCRD=PTRCRD+INDCRD+VALCRD
  855.  
  856. CALL OPENRB(NOMFIC,'mtx',CTITR,
  857. & TOTCRD,PTRCRD,INDCRD,VALCRD,
  858. & M,NVEC,NAUXD,0,
  859. & 'rua',PTRFMT,INDFMT,VALFMT)
  860. IF (IERR.NE.0) RETURN
  861. *
  862. * Conversion du stockage CSR (Compressed Sparse Row) vers
  863. * CSC (Compressed Sparse Column) = transposition
  864. * *******************************************************
  865. NTT=NBPTOT
  866. NJA=NVATOT
  867. NBVA=NJA
  868. SEGINI,KMOR2,KIZA2
  869. * CALL TRPMOR(NTT,NJA,KMORS.JA,KMORS.IA,KMOR2.JA,KMOR2.IA,
  870. * & 0,IRET)
  871. CALL TRMORS(NTT,NJA,KIZA.A,KMORS.JA,KMORS.IA,
  872. $ KIZA2.A,KMOR2.JA,KMOR2.IA)
  873. * Ecriture de la matrice
  874. WRITE(IOPER,PTRFMT) (KMOR2.IA(K),K=1,NBPTOT+1)
  875. WRITE(IOPER,INDFMT) (KMOR2.JA(K),K=1,NVATOT)
  876. WRITE(IOPER,VALFMT) (KIZA2.A(K),K=1,NVATOT)
  877. SEGSUP,KMOR2,KIZA2
  878. ENDIF
  879.  
  880. SEGSUP,KMORS,KIZA
  881.  
  882. ENDIF
  883.  
  884. CLOSE(UNIT=IOPER)
  885.  
  886.  
  887.  
  888. * +---------------------------------------------------------------+
  889. * | |
  890. * | G E O M E T R I E = F I C H I E R . G E O M |
  891. * | |
  892. * +---------------------------------------------------------------+
  893.  
  894. 200 CONTINUE
  895. IF (.NOT.ZGEOM) THEN
  896. OPEN(UNIT=IOPER,FILE=NOMFIC(1:LONG(NOMFIC))//'.geom'//CEXTN)
  897. CLOSE(UNIT=IOPER,STATUS='DELETE')
  898. GOTO 300
  899. ENDIF
  900.  
  901. IF (ITYP.EQ.2) THEN
  902. MNOEUD=IGEOMA
  903. SEGACT,MNOEUD
  904. NPOTOT=MNOEUD.NUM(/2)
  905. ENDIF
  906.  
  907. MYFMT='(E25.16)'
  908.  
  909. * Création du fichier .geom.mm
  910. * ****************************
  911. IF (IEXTN.EQ.1) THEN
  912. M=NPOTOT
  913. N=IDIM1
  914. NNZER=0
  915. CALL OPENMM(NOMFIC,'geom',CTITR,
  916. & M,N,NNZER,'matrix array real general')
  917.  
  918. * Création du fichier .geom.rb
  919. * ****************************
  920. ELSEIF (IEXTN.EQ.2) THEN
  921. M=NPOTOT
  922. NVEC=IDIM1
  923. NAUXD=0
  924. CALL OPENRB(NOMFIC,'geom',CTITR,
  925. & M,NVEC,NAUXD,0,
  926. & 0,0,0,0,
  927. & 'geos r',MYFMT,' ',' ')
  928. IF (IERR.NE.0) RETURN
  929. ENDIF
  930.  
  931. * Ecriture des coordonnées du maillage
  932. WRITE(IOPER,FMT=MYFMT)
  933. &((XCOOR((MNOEUD.NUM(1,K)-1)*IDIM1+J),K=1,NPOTOT),J=1,IDIM)
  934. WRITE(IOPER,FMT='(I12)')
  935. &(MNOEUD.NUM(1,K),K=1,NPOTOT)
  936.  
  937. IF (ITYP.EQ.2) SEGDES,MNOEUD
  938. CLOSE(UNIT=IOPER)
  939.  
  940.  
  941.  
  942. * +---------------------------------------------------------------+
  943. * | |
  944. * | I N C O N N U E S = F I C H I E R . I N C O |
  945. * | |
  946. * +---------------------------------------------------------------+
  947.  
  948. 300 CONTINUE
  949. IF (.NOT.ZINCO) THEN
  950. OPEN(UNIT=IOPER,FILE=NOMFIC(1:LONG(NOMFIC))//'.inco'//CEXTN)
  951. CLOSE(UNIT=IOPER,STATUS='DELETE')
  952. GOTO 400
  953. ENDIF
  954.  
  955. IF (ITYP.EQ.2) THEN
  956. IDDLPR=IINCPO
  957. IDDLDU=IDUAPO
  958. SEGACT,IDDLPR,IDDLDU
  959.  
  960. MIMIK=IIMIK
  961. MIDUA=IIDUA
  962. SEGACT,MIMIK,MIDUA
  963. ENDIF
  964.  
  965. NBVAR=NBPTOT+NBDTOT
  966. NBWR=2*NBVAR
  967. SEGINI,MIDATA
  968.  
  969. NBCOPR=IDDLPR.INCPO(/1)
  970. NBNOPR=IDDLPR.INCPO(/2)
  971. NBCODU=IDDLDU.INCPO(/1)
  972. NBNODU=IDDLDU.INCPO(/2)
  973.  
  974. * Inconnues primales
  975. DO J=1,NBCOPR
  976. DO I=1,NBNOPR
  977. IINC=IDDLPR.INCPO(J,I)
  978. IF (IINC.GT.0) THEN
  979. IWRIT(IINC)=I
  980. IWRIT(NBVAR+IINC)=J
  981. ENDIF
  982. ENDDO
  983. ENDDO
  984.  
  985. * Inconnues duales
  986. DO J=1,NBCODU
  987. DO I=1,NBNODU
  988. IINC=IDDLDU.INCPO(J,I)
  989. IF (IINC.GT.0) THEN
  990. IWRIT(NBPTOT+IINC)=I
  991. IWRIT(NBVAR+NBPTOT+IINC)=J
  992. ENDIF
  993. ENDDO
  994. ENDDO
  995.  
  996.  
  997. * Création du fichier .inco.mm
  998. * ****************************
  999. IF (IEXTN.EQ.1) THEN
  1000. M=NBVAR
  1001. N=2
  1002. NNZER=0
  1003. MYFMT='(I12)'
  1004. CALL OPENMM(NOMFIC,'inco',CTITR,
  1005. & M,N,NNZER,'matrix array real general')
  1006.  
  1007. * On ajoute le nom des composantes dans l'entête...
  1008. BACKSPACE(UNIT=IOPER)
  1009. WRITE(UNIT=IOPER,FMT='("%",/,A)')
  1010. & '% COMPOSANTES'
  1011. WRITE(UNIT=IOPER,FMT='("% PRIM ",I6," ",A4)')
  1012. & (J,IMIK(J),J=1,NBCOPR)
  1013. WRITE(UNIT=IOPER,FMT='("% DUAL ",I6," ",A4)')
  1014. & (J,IDUA(J),J=1,NBCODU)
  1015.  
  1016. * ...puis on reecrit ce que l'on avait effacé
  1017. WRITE(UNIT=IOPER,FMT='("%")')
  1018. WRITE(IOPER,FMT='(I12,1X,I12)') M,N
  1019.  
  1020. * Ecriture des noeud et composante associés à chaque variable
  1021. WRITE(IOPER,FMT=MYFMT) (IWRIT(K),K=1,NBWR)
  1022.  
  1023.  
  1024. * Création du fichier .inco.rb
  1025. * ****************************
  1026. ELSEIF (IEXTN.EQ.2) THEN
  1027. M=NBVAR
  1028. NVEC=2
  1029. NAUXD=0
  1030. MYFMT='(I12)'
  1031. CALL OPENRB(NOMFIC,'inco',CTITR,
  1032. & M,NVEC,NAUXD,0,
  1033. & 0,0,0,0,
  1034. & 'avl r',MYFMT,' ',' ')
  1035. IF (IERR.NE.0) RETURN
  1036.  
  1037. * Ecriture des noeud et composante associés à chaque variable
  1038. WRITE(IOPER,FMT=MYFMT) (IWRIT(K),K=1,NBWR)
  1039.  
  1040. * On ajoute le nom des composantes en fin de fichier
  1041. WRITE(UNIT=IOPER,FMT='("%",/,A)')
  1042. & '% COMPOSANTES'
  1043. WRITE(UNIT=IOPER,FMT='("% PRIM ",I6," ",A4)')
  1044. & (J,IMIK(J),J=1,NBCOPR)
  1045. WRITE(UNIT=IOPER,FMT='("% DUAL ",I6," ",A4)')
  1046. & (J,IDUA(J),J=1,NBCODU)
  1047.  
  1048. ENDIF
  1049.  
  1050.  
  1051. CLOSE(UNIT=IOPER)
  1052. IF (ITYP.EQ.2) SEGDES,MIMIK,MIDUA,IDDLPR,IDDLDU
  1053. SEGSUP,MIDATA
  1054.  
  1055.  
  1056.  
  1057.  
  1058.  
  1059.  
  1060. * +---------------------------------------------------------------+
  1061. * | |
  1062. * | S E C O N D - M E M B R E = F I C H I E R . R H S |
  1063. * | |
  1064. * +---------------------------------------------------------------+
  1065.  
  1066. 400 CONTINUE
  1067. IF (.NOT.ZFORC) THEN
  1068. OPEN(UNIT=IOPER,FILE=NOMFIC(1:LONG(NOMFIC))//'.rhs'//CEXTN)
  1069. CLOSE(UNIT=IOPER,STATUS='DELETE')
  1070. GOTO 500
  1071. ENDIF
  1072.  
  1073. MYFMT='(E25.16)'
  1074.  
  1075. * Création du fichier .rhs.mm
  1076. * ****************************
  1077. IF (IEXTN.EQ.1.AND.ISMBR.NE.0) THEN
  1078. M=NBDTOT
  1079. N=1
  1080. NNZER=0
  1081. CALL OPENMM(NOMFIC,'rhs',CTITR,
  1082. & M,N,NNZER,'matrix array real general')
  1083.  
  1084. * Ecriture des valeur du second membre
  1085. WRITE(IOPER,FMT=MYFMT)
  1086. & (ISMBR.VECTBB(K),K=1,NBDTOT)
  1087. SEGSUP,ISMBR
  1088. ELSE
  1089. WRITE(IOIMP,*) 'Pas de sortie FORC pour le moment'
  1090. CALL ERREUR(251)
  1091. ENDIF
  1092.  
  1093. IF (IERR.NE.0) RETURN
  1094.  
  1095.  
  1096. * +---------------------------------------------------------------+
  1097. * | |
  1098. * | C O N N E C T I V I T É = F I C H I E R . C O N N |
  1099. * | |
  1100. * +---------------------------------------------------------------+
  1101.  
  1102. 500 CONTINUE
  1103. IF (.NOT.ZELEM) GOTO 600
  1104. WRITE(IOIMP,*) 'Pas de sortie CONN pour le moment'
  1105. CALL ERREUR(251)
  1106. IF (IERR.NE.0) RETURN
  1107.  
  1108.  
  1109.  
  1110. * +---------------------------------------------------------------+
  1111. * | |
  1112. * | S O L U T I O N D E R E F E R E N C E |
  1113. * | = F I C H I E R . S O L U |
  1114. * | |
  1115. * +---------------------------------------------------------------+
  1116.  
  1117. 600 CONTINUE
  1118. IF (.NOT.ZSOLU) GOTO 700
  1119. WRITE(IOIMP,*) 'Pas de sortie SOLU pour le moment'
  1120. CALL ERREUR(251)
  1121. IF (IERR.NE.0) RETURN
  1122.  
  1123.  
  1124. * +---------------------------------------------------------------+
  1125. * | |
  1126. * | R É S U L T A T D E C A L C U L |
  1127. * | = F I C H I E R . R E S U |
  1128. * | |
  1129. * +---------------------------------------------------------------+
  1130.  
  1131. 700 CONTINUE
  1132. IF (.NOT.ZRESU) THEN
  1133. OPEN(UNIT=IOPER,FILE=NOMFIC(1:LONG(NOMFIC))//'.resu'//CEXTN)
  1134. CLOSE(UNIT=IOPER,STATUS='DELETE')
  1135. GOTO 999
  1136. ENDIF
  1137.  
  1138. MYFMT='(E25.16)'
  1139.  
  1140. * Création du fichier .resu.mm
  1141. * ****************************
  1142. IF (IEXTN.EQ.1.AND.INCX.NE.0) THEN
  1143. M=NBPTOT
  1144. N=1
  1145. NNZER=0
  1146. CALL OPENMM(NOMFIC,'resu',CTITR,
  1147. & M,N,NNZER,'matrix array real general')
  1148.  
  1149. * Ecriture des valeur du second membre
  1150. WRITE(IOPER,FMT=MYFMT)
  1151. & (INCX.VECTBB(K),K=1,NBPTOT)
  1152. SEGSUP,INCX
  1153. ELSE
  1154. WRITE(IOIMP,*) 'Pas de sortie RESU pour le moment'
  1155. CALL ERREUR(251)
  1156. ENDIF
  1157. IF (IERR.NE.0) RETURN
  1158.  
  1159.  
  1160. * +---------------------------------------------------------------+
  1161. * | |
  1162. * | M É N A G E E T S O R T I E |
  1163. * | |
  1164. * +---------------------------------------------------------------+
  1165.  
  1166. 999 CONTINUE
  1167.  
  1168.  
  1169. IF (ITYP.EQ.1) SEGSUP,MNOEUD,MIMIK,MIDUA,IDDLPR,IDDLDU
  1170. IF (ITYP.EQ.2) SEGDES,MMATRI
  1171. SEGDES,MRIGID
  1172. RETURN
  1173.  
  1174. 666 CONTINUE
  1175. WRITE(IOIMP,*) 'Lagrangian multiplier detected. Untreated case'
  1176. MOTERR(1:8)='SORMAT '
  1177. CALL ERREUR(349)
  1178. RETURN
  1179.  
  1180. END
  1181.  
  1182.  
  1183.  
  1184.  
  1185.  
  1186.  
  1187.  
  1188.  
  1189.  
  1190.  

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