Télécharger la plaquette Cast3M


fr en

Télécharger trialu.eso

Retour à la liste

Numérotation des lignes : 

trialu

  1. C TRIALU    SOURCE    GOUNAND   22/08/25    21:15:12     11434          
  2.       SUBROUTINE TRIALU(MATRIK,AMORS,AISA,
  3.      $     IDMAT,
  4.      $     IMPR,IRET)
  5.       IMPLICIT REAL*8 (A-H,O-Z)
  6.       IMPLICIT INTEGER (I-N)
  7. C************************************************************************
  8. C
  9. C    Triangulation d'une matrice symetrique ou non
  10. C    On decoupe la matrice en bloc
  11. C    Si IDMAT ne 0 la matrice est deja triangulee on saute
  12. C
  13. C HISTORIQUE : 22/03/00 gounand
  14. C Amélioration de la détection des pivots nulles (pas encore parfaite)
  15. C HISTORIQUE : 13/01/00 Adaptation au nouvel assemblage :
  16. C      Modifications cosmétiques.
  17. C      Correction sur la détection des diagonales nulles (pas encore
  18. C      parfaite).
  19. C      Modification de la taille des paquets pour tenir dans le cache
  20. C      (dépendant machine).
  21. C      Inlining manuel des produits scalaires de vecteurs (sdot) qui
  22. C      sont appelés de nombreuses fois...
  23. C      Gestion correcte des allocations mémoire.
  24. C HISTORIQUE :
  25. C HISTORIQUE :
  26. C************************************************************************
  27. *      REAL*8 SDOT
  28.  
  29. -INC PPARAM
  30. -INC CCOPTIO
  31. -INC CCREEL
  33.       INTEGER NTT,NPT,NBLK,KTRING
  34.       INTEGER NBVA
  35.       POINTEUR ISAU.IZA,ISAL.IZA
  36.       POINTEUR ISAU0.IZA,ISAL0.IZA
  37.       POINTEUR ISAU00.IZA,ISAL00.IZA
  38.       POINTEUR AMORS.PMORS
  39.       POINTEUR AISA.IZA
  40.       POINTEUR PMTRAN.PMORS
  41.       POINTEUR MYMINC.MINC
  42. -INC SMELEME
  43.       POINTEUR MYPTS.MELEME
  44. -INC SMLREEL
  45.       POINTEUR GGIS.MLREEL
  46. *STAT-INC SMSTAT
  47.       INTEGER IMPR,IRET
  48. *
  49.       INTEGER NBVD
  50.       SEGMENT IZD
  51.       REAL*8 D(NBVD)
  52.       ENDSEGMENT
  53.       INTEGER IBLSIZ
  54. *      REAL*8 GGIS
  55.       LOGICAL LACTL,LACTU
  56. * Le logique LISDDL indique si l'on est capable de faire
  57. * la correspondance numéro DDL <-> numéro de point, nom d'inconnue
  58.       LOGICAL LISDDL
  59.       INTEGER TRAMIN
  60.       INTEGER I,ISTRT,ISTOP
  61.       INTEGER J,JSTRT,JSTOP
  62.       INTEGER NTTDDL
  63.       REAL*8 TNORM,YZPREC,US
  64. C*******************************************************************
  65. C
  66. C
  67.       IF (IMPR.GT.5) WRITE(IOIMP,*) 'Entrée dans trialu'
  68.       IKOMPT=0
  69. *      GGIS=XZPREC
  70.       SEGACT MATRIK
  71.       KTRING=MATRIK.KIDMAT(9)
  72.       KSIM=MATRIK.KSYM
  73.       IDMAT=MATRIK.KIDMAT(1)
  74.       MYMINC=MATRIK.KMINC
  75.       MYPTS=MATRIK.KISPGT
  76.       SEGDES MATRIK
  77.       IF(KTRING.NE.2) THEN
  78. *STAT         CALL INMSTA(MSTAT,1)
  79.          SEGACT MATRIK*MOD
  80.          MATRIK.KIDMAT(9)=2
  81.          SEGDES MATRIK
  82.          SEGACT IDMAT*MOD
  83. C Initialisation des KZA
  84.          SEGACT AMORS
  85.          NTT=AMORS.IA(/1)-1
  86.          NPT=IDMAT.NUAN(/1)
  87.          LISDDL=(MYMINC.NE.0.AND.MYPTS.NE.0.AND.NPT.GT.0)
  88.          NBLK=0
  89.          SEGADJ,IDMAT
  90.          DO 9 ITT=1,NTT
  91.             IDMAT.KZA(ITT)=ITT-AMORS.JA(AMORS.IA(ITT))+1
  92.  9       CONTINUE
  93.          SEGDES AMORS
  94. C On construit la transposée du profil
  95. C (La matrice Skyline est supposée avoir une structure symétrique)
  96. C La transposée du profil n'a pas ses colonnes ordonnées...
  97.          CALL MAKPMT(AMORS,
  98.      $        PMTRAN,
  99.      $        IMPR,IRET)
  100.          IF (IRET.NE.0) GOTO 9999
  101.          SEGACT PMTRAN
  102.          DO 91 ITT=1,NTT
  103.             JSTRT=PMTRAN.IA(ITT)
  104.             JSTOP=PMTRAN.IA(ITT+1)-1
  105.             TRAMIN=PMTRAN.JA(JSTRT)
  106.             DO 912 J=JSTRT+1,JSTOP
  107.                TRAMIN=MIN(TRAMIN,PMTRAN.JA(J))
  108.  912        CONTINUE
  109.             TRAMIN=DIM(ITT,TRAMIN)+1
  110.             IDMAT.KZA(ITT)=MAX(IDMAT.KZA(ITT),TRAMIN)
  111.  91      CONTINUE
  112.          SEGSUP PMTRAN
  113. *STAT         CALL PRMSTA('Init. KZA',MSTAT,1)
  114. C
  115. C On cree les blocs
  116. C
  117. C La matrice est stockée par blocs de plusieurs lignes. En effet, il faut
  118. C l'équivalent de 900 opérations pour activer/désactiver un segment. Si
  119. C une ligne de la matrice comprend 50 éléments, on fait 100 opérations
  120. C utiles et 900 inutiles lors de la résolution.
  121. C La taille des blocs est fixée par le paramètre IBLSIZ. La diagonale est
  122. C stockée dans un segment à part.
  123. C
  124. C Il nous faut constituer un certain nombre de tableaux. Le premier segment
  125. C décrit la matrice. Il contient un tableau donnant les numéros des premières
  126. C lignes de chaque bloc et un tableau de pointeurs sur des segments
  127. C descripteurs de bloc.
  128. C
  129. C Calcul de la longueur des blocs, on choisit ce qu'il faut
  130. C
  131. C Attention !! il y a au moins deux blocs activés en même temps et au plus 4
  132. C Peut-être qu'une bonne formule fait intervenir la taille du cache...
  133. C     (Cache-Place(KZA,NUIA,DESCL...)) / 4 (blocs simultanés)
  134. *         IBLSIZ=OOOVAL(1,1)/100000
  135. *  Pour 128Ko de cache => 8 KMots
  136.          IBLSIZ=8000
  137. C
  138. C Il faut d'abord déterminer le nombre de blocs.
  139. C
  140.          LONG=0
  141.          NTT=KZA(/1)
  142.          DO 96 I=1,NTT
  143.             LONG=LONG+KZA(I)
  144.  96      CONTINUE
  145.          NBLK=1
  146.          IF(LONG-NTT.GT.IBLSIZ) THEN
  147.             ILBL=0
  148.             DO 12 KK=2,NTT
  149.                LA=KZA(KK)-1
  150.                IF(LA.GT.IBLSIZ) THEN
  151.                   NBLK=NBLK+2
  152.                   ILBL=0
  153.                ELSE
  154.                   ILBL=ILBL+LA
  155.                   IF(ILBL.GT.IBLSIZ) THEN
  156.                      NBLK=NBLK+1
  157.                      ILBL=LA
  158.                   ENDIF
  159.                ENDIF
  160.  12         CONTINUE
  161.          ENDIF
  162. C
  163. C On peut allouer le segment descripteur de la matrice. On va pouvoir
  164. C remplir des numéros de ligne des débuts de bloc.
  165. C
  166.          NTT=KZA(/1)
  167.          NPT=NUAN(/1)
  168.          SEGADJ IDMAT
  169.          IF(NBLK.EQ.1) THEN
  170.             NLDBLK(1)=1
  171.             NLDBLK(2)=NTT+1
  172.          ELSE
  173.             IBLK=1
  174.             ILBL=0
  175.             NLDBLK(1)=1
  176.             NLDBLK(NBLK+1)=NTT+1
  177.             DO 13 KK=2,NTT
  178.                LA=KZA(KK)-1
  179.                IF(LA.GT.IBLSIZ) THEN
  180.                   NLDBLK(IBLK+1)=KK
  181.                   NLDBLK(IBLK+2)=KK+1
  182.                   IBLK=IBLK+2
  183.                   ILBL=0
  184.                ELSE
  185.                   ILBL=ILBL+LA
  186.                   IF(ILBL.GT.IBLSIZ) THEN
  187.                      NLDBLK(IBLK+1)=KK
  188.                      ILBL=LA
  189.                      IBLK=IBLK+1
  190.                   ENDIF
  191.                ENDIF
  192.  13         CONTINUE
  193.          ENDIF
  194.          IF(IMPR.GE.2) THEN
  195.             IF(KSIM.EQ.0) THEN
  196.                WRITE(IOIMP,*) '* NOMBRE DE BLOCS DE ',IBLSIZ,
  197.      $              ' MOTS   : ',NBLK
  198.             ELSE
  199.                WRITE(IOIMP,*)'* NOMBRE DE BLOCS DE ',IBLSIZ,
  200.      $              ' MOTS   : 2 X ',NBLK
  201.             ENDIF
  202.          ENDIF
  203. C
  204. C On Remplit les blocs
  205. C
  206.          SEGACT AMORS
  207.          SEGACT AISA
  208. *  On construit le vecteurs des normes des lignes
  209.          NTTDDL=AMORS.IA(/1)-1
  210.          JG=NTTDDL
  211.          SEGINI GGIS
  212.          YZPREC=XZPREC**(0.9D0)
  213. *         YZPREC=1.D-9
  214.          DO 152 ITTDDL=1,NTTDDL
  215. * Calcul de la norme de la ittème ligne
  216.             ISTRT=AMORS.IA(ITTDDL)
  217.             ISTOP=AMORS.IA(ITTDDL+1)-1
  218.             TNORM=XZERO
  219.             DO 1521 IJA=ISTRT,ISTOP
  220.                TNORM=TNORM+ABS(AISA.A(IJA))
  221.  1521       CONTINUE
  222. *            WRITE(IOIMP,*) 'ITTDDL=',ITTDDL,' TNORM=',TNORM
  223. *
  224. C            IF (TNORM.EQ.0.D0) THEN
  225. C               WRITE(IOIMP,*) 'The ',ITTDDL,'th line is zero'
  226. CC Résolution impossible détectée au noeud %i1 pour l'inconnue %m1:4
  227. C               SEGACT MYMINC
  228. C               SEGACT MYPTS
  229. C               N2=IDMAT.NUNA(ITTDDL)
  230. C               CALL DDL2PI(N2,MYMINC,
  231. C     $              IPT,IBI,
  232. C     $              IMPR,IRET)
  233. C               IF (IRET.NE.0) GOTO 9999
  234. C               INTERR(1)=MYPTS.NUM(1,IPT)
  235. C               MOTERR(1:8)=MYMINC.LISINC(IBI)
  236. C               CALL ERREUR(143)
  237. C               SEGDES MYPTS
  238. C               SEGDES MYMINC
  239. C               GOTO 9999
  240. C            ENDIF
  241.             GGIS.PROG(ITTDDL)=(TNORM/DBLE(ISTOP-ISTRT+1))*YZPREC
  242.  152     CONTINUE
  243.          NBVD=KZA(/1)
  244.          SEGINI IZD
  245.          IDMAT.IDIAG=IZD
  246.          ISAU=0
  247.          DO 601 IBLK=1,NBLK
  248. * Gestion du CTRL-C
  249.             if (ierr.NE.0) return
  250. C NUMÉROS DES LIGNES DE DÉBUT ET DE FIN DE BLOC KJD,KJF
  251.             KJD=NLDBLK(IBLK)
  252.             KJF=NLDBLK(IBLK+1)-1
  253.             KA=0
  254.             KA=KA+KZA(1)-1
  255.             NUIA(1,1)=1
  256.             DO 600 KJ=KJD,KJF
  257.                NUIA(KJ,1)=IBLK
  258.                NUIA(KJ,2)=KA
  259.                KA=KA+KZA(KJ)-1
  260.  600        CONTINUE
  261. C ALLOCATION DE LA MÉMOIRE POUR LE BLOC
  262.             NBVA=KA
  263. *            WRITE(IOIMP,*) '* BLOC LOWER NUMERO : ',IBLK,' TAILLE = '
  264. *     $           ,NBVA,' MOTS'
  265.             SEGINI IZA
  266.             IDESCL(IBLK)=IZA
  267.             IDESCU(IBLK)=IZA
  268.             IF(KSIM.EQ.2)IDESCU(IBLK)=NBVA
  269.             IF(IBLK.EQ.1)THEN
  270.                KJD=1
  271.                D(1)=AISA.A(1)
  272.             ENDIF
  273.             DO 602 LI=KJD,KJF
  274.                NBK=AMORS.IA(LI+1)-AMORS.IA(LI)
  275.                I1=AMORS.IA(LI)
  276.                DO 603 K=1,NBK
  277.                   MA=I1+K-1
  278.                   KO=AMORS.JA(MA)
  279.                   IF(LI.GT.KO)THEN
  280.                      NA=KO-LI+KZA(LI)+NUIA(LI,2)
  281.                      A(NA)=AISA.A(MA)
  282.                   ELSEIF(LI.EQ.KO.AND.LI.NE.1)THEN
  283.                      D(LI)=AISA.A(MA)
  284.                   ELSE
  285.                      GO TO 602
  286.                   ENDIF
  287.  603           CONTINUE
  288.  602        CONTINUE
  289.             SEGDES IZA
  290.  601     CONTINUE
  291.          IF(KSIM.EQ.2)THEN
  292.             DO 604 IBLK=1,NBLK
  293. * Gestion du CTRL-C
  294.                if (ierr.NE.0) return
  295.                KJD=NLDBLK(IBLK)
  296.                KJF=NLDBLK(IBLK+1)-1
  297.                NBVA=IDESCU(IBLK)
  298.                SEGINI IZA
  299.                IDESCU(IBLK)=IZA
  300.                IF(IBLK.EQ.1)THEN
  301.                   KJD=1
  302.                ENDIF
  303.                LDEB=KJD
  304.                DO 605 KO=KJD,KJF
  305.                   L0=KO-KZA(KO)+1
  306.                   IF(LDEB.GT.L0)LDEB=L0
  307.  605           CONTINUE
  308.                DO 606 LI=LDEB,KJF
  309.                   NBK=AMORS.IA(LI+1)-AMORS.IA(LI)
  310.                   I1=AMORS.IA(LI)
  311.                   DO 607 K=NBK,1,-1
  312.                      MA=I1+K-1
  313.                      KO=AMORS.JA(MA)
  314.                      IF(KO.GT.KJF)GO TO 607
  315.                      IF(KO.LT.KJD)GO TO 606
  316.                      IF(LI.LT.KO)THEN
  317.                         NA=LI-KO+KZA(KO)+NUIA(KO,2)
  318.                         A(NA)=AISA.A(MA)
  319.                      ELSE
  320.                         GO TO 606
  321.                      ENDIF
  322.  607              CONTINUE
  323.  606           CONTINUE
  324.                SEGDES IZA
  325.  604        CONTINUE
  326.          ENDIF
  327.          SEGDES AISA
  328.          SEGDES AMORS
  329. *STAT         CALL PRMSTA('Remplissage des blocs',MSTAT,1)
  330.          IF(KSIM.EQ.0.OR.KSIM.EQ.2)THEN
  331.             DO 200 IBLK=1,NBLK
  332. * Gestion du CTRL-C
  333.                if (ierr.NE.0) return
  334.                KJD=NLDBLK(IBLK)
  335.                KJF=NLDBLK(IBLK+1)-1
  336.                ISAU=IDESCU(IBLK)
  337.                ISAL=IDESCL(IBLK)
  338.                SEGACT ISAU*MOD
  339.                SEGACT ISAL*MOD
  340.                IF(IBLK.NE.1)THEN
  341.                   ISAL0=0
  342.                   ISAU0=0
  343. C     SEGACT ISAL0
  344. C     SEGACT ISAU0
  345.                   ISAL00=0
  346.                   ISAU00=0
  347.                   LACTL=.FALSE.
  348.                   LACTU=.FALSE.
  349.                   KJD0=KJD-KZA(KJD)+1
  350.                   DO 304 I=KJD+1,KJF
  351.                      KLD=I-KZA(I)+1
  352.                      IF(KLD.LT.KJD0)KJD0=KLD
  353.  304              CONTINUE
  354.                   DO 302 I=KJD0,KJD-1
  355.                      LDEB=I-KZA(I)+1
  356.                      DO 301 N=KJD,KJF
  357.                         KDEB=N-KZA(N)+1
  358.                         IF(KDEB.LE.I) THEN
  359. C Calcul colonne n
  360.                            K1=MAX(KDEB,LDEB)
  361.                            LG=I-K1
  362.                            KL=NUIA(I,2)+K1-LDEB+1
  363.                            KU=NUIA(N,2)+K1-KDEB+1
  364.                            ISAL00=IDESCL(NUIA(I,1))
  365.                            LACTL=(ISAL00.NE.ISAL0.AND.ISAL00.NE.ISAL)
  366.                            IF (LACTL) THEN
  367.                               IF (ISAL0.NE.0) THEN
  368.                                  SEGDES ISAL0
  369.                               ENDIF
  370.                               SEGACT ISAL00
  371.                            ENDIF
  372.                            ISAL0=ISAL00
  373. *                           US=SDOT(LG,ISAL00.A(KL),1,ISAU.A(KU),1)
  374.                            US=0.D0
  375.                            DO ILG=0,LG-1
  376.                               US=US+ISAL00.A(KL+ILG)*ISAU.A(KU+ILG)
  377.                            ENDDO
  378.                            KU=NUIA(N,2)+I-KDEB+1
  379.                            ISAU.A(KU)=(ISAU.A(KU)-US)
  380. C     WRITE(IOIMP,*)' KU=',KU,' I=',I,' A(KU)=',ISAU.A(KU)
  381.                            IF(KSIM.NE.0)THEN
  382. C Calcul ligne n
  383.                               KU=NUIA(I,2)+K1-LDEB+1
  384.                               KL=NUIA(N,2)+K1-KDEB+1
  385.                               ISAU00=IDESCU(NUIA(I,1))
  386.                               LACTU=(ISAU00.NE.ISAU0.AND.ISAU00.NE.ISAU)
  387.                               IF (LACTU) THEN
  388.                                  IF (ISAU0.NE.0) THEN
  389.                                     SEGDES ISAU0
  390.                                  ENDIF
  391.                                  SEGACT ISAU00
  392.                               ENDIF
  393.                               ISAU0=ISAU00
  394. *                              US=SDOT(LG,ISAL.A(KL),1,ISAU00.A(KU),1)
  395.                               US=0.D0
  396.                               DO ILG=0,LG-1
  397.                                  US=US+ISAL.A(KL+ILG)*ISAU00.A(KU+ILG)
  398.                               ENDDO
  399.                               KL=NUIA(N,2)+I-KDEB+1
  400.                               ISAL.A(KL)=(ISAL.A(KL)-US)
  401.                            ENDIF
  402.                         ENDIF
  403.  301                 CONTINUE
  404.  302              CONTINUE
  405.                   IF (ISAL00.NE.ISAL.AND.ISAL00.NE.0) THEN
  406.                      SEGDES ISAL00
  407.                   ENDIF
  408.                   IF (ISAU00.NE.ISAU.AND.ISAU00.NE.0) THEN
  409.                      SEGDES ISAU00
  410.                   ENDIF
  411.                ENDIF
  412.                DO 202 I=KJD,KJF
  413.                   LDEB=I-KZA(I)+1
  414.                   KL=NUIA(I,2)
  415.                   DO 2011 II=LDEB,I-1
  416.                      ISAL.A(KL+II-LDEB+1)=ISAL.A(KL+II-LDEB+1)/D(II)
  417.  2011             CONTINUE
  418.                   IF(KSIM.EQ.2)THEN
  419.                      KL=NUIA(I,2)
  420.                      DO 2012 II=LDEB,I-1
  421.                         ISAU.A(KL+II-LDEB+1)=ISAU.A(KL+II-LDEB+1)/D(II)
  422.  2012                CONTINUE
  423.                   ENDIF
  424. C Calcul diagonale
  425.                   N=I
  426.                   IF(N.NE.1) THEN
  427.                      LG=N-LDEB
  428.                      KK=NUIA(N,2)+1
  429.                      US=0.D0
  430.                      DO 3003 II=1,LG
  431.                         US=US+
  432.      $                       ISAL.A(KK+II-1)
  433.      $                       *D(II+LDEB-1)*ISAU.A(KK+II-1)
  434.  3003                CONTINUE
  435.                      D(N)=D(N)-US
  436.                   ENDIF
  437. C! DEBUT
  438.                   ADK=ABS(D(N))
  439. *                  WRITE(IOIMP,*) 'N=',N,' ADK=',ADK
  440.                   IF(ADK.LE.GGIS.PROG(N))THEN
  441.                      IKOMPT=IKOMPT+1
  442. *
  443. * Désactivation temporaire, sinon des cas-tests (mal conçus !) plantent
  444. *
  445. C Résolution impossible détectée au noeud %i1 pour l'inconnue %m1:4
  446. *                     INTERR(1)=N
  447. *                     MOTERR(1: 8) = '        '
  448. *                     CALL ERREUR(143)
  449.                      IF (IKOMPT.EQ.1) THEN
  450.                         WRITE(IOIMP,*) '!!!!!!!!!!!!! WARNING !!!!!!!!!'
  451.                         IF (LISDDL) THEN
  452.                            SEGACT MYMINC
  453.                            SEGACT MYPTS
  454.                            DO IDDL=1,MIN(10,NTTDDL)
  455.                               N2=IDMAT.NUNA(IDDL)
  456.                               CALL DDL2PI(N2,MYMINC,
  457.      $                             IPT,IBI,
  458.      $                             IMPR,IRET)
  459.                               IF (IRET.NE.0) GOTO 9999
  460.                               WRITE(IOIMP,*) 'ddl ',IDDL
  461.      $                             ,' = point numero ',MYPTS.NUM(1,IPT),
  462.      $                             ' inconnue ',MYMINC.LISINC(IBI)
  463.                            ENDDO
  464.                            SEGDES MYPTS
  465.                            SEGDES MYMINC
  466.                         ENDIF
  467.                      ENDIF
  468.                      IF (IKOMPT.LE.10.OR.ADK.EQ.0.D0) THEN
  469.                         WRITE(IOIMP,*) 'Au ddl ',N,'le pivot vaut ',ADK,
  470.      $                       ' ; norme de la ligne =',GGIS.PROG(N)
  471.      $                       /YZPREC
  472.                         IF (LISDDL) THEN
  473.                            N2=IDMAT.NUNA(N)
  474.                            CALL DDL2PI(N2,MYMINC,
  475.      $                          IPT,IBI,
  476.      $                          IMPR,IRET)
  477.                            IF (IRET.NE.0) GOTO 9999
  478.                            SEGACT MYMINC
  479.                            SEGACT MYPTS
  480.                            IF (ADK.EQ.0.D0) THEN
  481.                               INTERR(1)=MYPTS.NUM(1,IPT)
  482.                               MOTERR(1:8)=MYMINC.LISINC(IBI)
  483.                               CALL ERREUR(143)
  484.                               GOTO 9999
  485.                            ELSE
  486.                               WRITE(IOIMP,*) 'ddl ',N,' = point numero '
  487.      $                             ,MYPTS.NUM(1,IPT),' inconnue ',MYMINC
  488.      $                             .LISINC(IBI)
  489.                            ENDIF
  490.                            SEGDES MYPTS
  491.                            SEGDES MYMINC
  492.                         ELSE
  493.                            IF (ADK.EQ.0.D0) THEN
  494.                               INTERR(1)=N
  495. *  49 2
  496. *Matrice singulière. Numéro de ligne =%i1
  497.                               CALL ERREUR(49)
  498.                               GOTO 9999
  499.                            ENDIF
  500.                         ENDIF
  501.                         IF (ADK.EQ.0.D0) RETURN
  502.                      ELSE IF (IKOMPT.EQ.11) THEN
  503.                         WRITE(IOIMP,*) 'Plus de 10 pivots petits'
  504.                      ELSE
  505. * Temporairement désactivé pour cause de mauvais cas-tests : difasyk2Dax.dgibi
  506. *                        GOTO 9999
  507.                      ENDIF
  508.                      D(N)=GGIS.PROG(N)
  509. *                     D(N)=1.D-200
  510.                   ENDIF
  511. C! FIN
  512.                   I1=I+1
  513.                   DO 201 N=I1,KJF
  514.                      KDEB=N-KZA(N)+1
  515.                      IF(KDEB.LE.I) THEN
  516. C Calcul colonne n
  517.                         K1=MAX(KDEB,LDEB)
  518.                         LG=I-K1
  519.                         KL=NUIA(I,2)+K1-LDEB+1
  520.                         KU=NUIA(N,2)+K1-KDEB+1
  521. *                        US=SDOT(LG,ISAL.A(KL),1,ISAU.A(KU),1)
  522.                         US=0.D0
  523.                         DO ILG=0,LG-1
  524.                            US=US+ISAL.A(KL+ILG)*ISAU.A(KU+ILG)
  525.                         ENDDO
  526. *
  527.                         KU=NUIA(N,2)+I-KDEB+1
  528.                         ISAU.A(KU)=(ISAU.A(KU)-US)
  529.                         IF(KSIM.NE.0) THEN
  530. C Calcul ligne n
  531.                            KU=NUIA(I,2)+K1-LDEB+1
  532.                            KL=NUIA(N,2)+K1-KDEB+1
  533. *                           US=SDOT(LG,ISAL.A(KL),1,ISAU.A(KU),1)
  534.                            US=0.D0
  535.                            DO ILG=0,LG-1
  536.                               US=US+ISAL.A(KL+ILG)*ISAU.A(KU+ILG)
  537.                            ENDDO
  538. *
  539.                            KL=NUIA(N,2)+I-KDEB+1
  540.                            ISAL.A(KL)=(ISAL.A(KL)-US)
  541.                         ENDIF
  542.                      ENDIF
  543.  201              CONTINUE
  544.  202           CONTINUE
  545.                SEGDES ISAL
  546.                SEGDES ISAU
  547.  200        CONTINUE
  548.             SEGSUP GGIS
  549.          ELSE
  550.             WRITE(IOIMP,*)' TRIALU : KSIM=',KSIM,' Cas non prevu '
  551.             GOTO 9999
  552.          ENDIF
  553. *STAT         CALL PRMSTA('Triangulation',MSTAT,1)
  554.          SEGDES IZD
  555.          SEGDES IDMAT
  556.       ENDIF
  557. *
  558. * Normal termination
  559. *
  560.       IRET=0
  561.       RETURN
  562. *
  563. * Format handling
  564. *
  565. *
  566. * Error handling
  567. *
  568.  9999 CONTINUE
  569.       IRET=1
  570.       WRITE(IOIMP,*) 'An error was detected in subroutine trialu.eso'
  571.       RETURN
  572. *
  573. * End of subroutine TRIALU
  574. *
  575.       END
  576.  
  577.  

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