1 : C EXDIAR SOURCE GOUNAND 12/11/12 21:15:03 7564 2 : SUBROUTINE EXDIAR(MRIGID,MCHPOI) 3 : IMPLICIT REAL*8 (A-H,O-Z) 4 : IMPLICIT INTEGER (I-N) 5 : C************************************************************************* 6 : C Operateur EXDIAR 7 : C 8 : C OBJET : Extrait la diagonale d'une matrice. 9 : C 10 : C*********************************************************************** 11 : C HISTORIQUE : 06/06/1 : Première version 12 : C 13 : C*********************************************************************** 14 : -INC CCOPTIO 15 : -INC CCHAMP 16 : -INC SMCHPOI 17 : -INC SMRIGID 18 : -INC SMELEME 19 : -INC SMCOORD 20 : -INC TMTRAV 21 : * 22 : segment trav1 23 : character*4 compd(nbincd),compdd(nbincd) 24 : endsegment 25 : segment trav2 26 : integer cor1p(nligrp),cor1d(nligrd) 27 : endsegment 28 : * 29 : * Pour chaque segment DESCR, on a un segment CORES2 où 30 : * CORDIA donne les couples (NLIGRP,NLIGRD) qui sont des 31 : * termes diagonaux avec les cas particuliers suivant : 32 : * CORES2=0 : il n'y a pas de termes diagonaux dans ces matrices 33 : * élémentaires 34 : * CORES2=-1 : les termes diagonaux de la matrices élémentaire (carrée) 35 : * sont des termes diagonaux de la matrice globale 36 : * on a aussi a segment CORES3 où COMTRA indique 37 : * le nom d'inconnue primal+numéro d'harmonique suivant 38 : * le repérage induit par MTRA1 et MTRA2 39 : * 40 : SEGMENT CORES1 41 : INTEGER ITRAV2(NRIGEL) 42 : INTEGER IPCOR2(NRIGEL) 43 : INTEGER IPCOR3(NRIGEL) 44 : ENDSEGMENT 45 : SEGMENT CORES2 46 : INTEGER CORDIA(2,NLIGP) 47 : ENDSEGMENT 48 : SEGMENT CORES3 49 : INTEGER COMTRA(NLIGP) 50 : ENDSEGMENT 51 : * 52 : SEGMENT ICPR(XCOOR(/1)/(IDIM+1)) 53 : SEGMENT MTRA1 54 : * CHARACTER*4 ICOMP(0) 55 : INTEGER ICOMP(0) 56 : ENDSEGMENT 57 : SEGMENT MTRA2 58 : INTEGER MHAR(0) 59 : ENDSEGMENT 60 : logical descar 61 : * 62 : * Construisons d'abord l'ensemble des noms d'inconnues 63 : * primales, puis les duales associées (segment trav1) 64 : * 65 : C WRITE(IOIMP,*) 'Coucou exdiar' 66 : SEGACT MRIGID 67 : nbincd=0 68 : segini trav1 69 : NRIGEL=COERIG(/1) 70 : DO IRIG=1,NRIGEL 71 : DESCR=IRIGEL(3,IRIG) 72 : SEGACT DESCR 73 : do 45 iligrp=1,lisinc(/2) 74 : do 55 inc=1,nbincd 75 : if (compd(inc).eq.lisinc(iligrp)) goto 45 76 : 55 continue 77 : nbincd=nbincd+1 78 : segadj trav1 79 : compd(nbincd)=lisinc(iligrp) 80 : 45 continue 81 : SEGDES DESCR 82 : ENDDO 83 : do iincd=1,nbincd 83 : CALL PLACE(NOMDD,LNOMDD,idx,compd(iincd)) 84 : IF (idx.EQ.0) THEN 85 : compdd(iincd)=compd(iincd) 86 : ELSE 87 : compdd(iincd)=NOMDU(idx) 88 : ENDIF 89 : ENDDO 90 : C WRITE(IOIMP,*) ' 1) compd' 91 : C WRITE (IOIMP,2019) (compd(I),I=1,compd(/2)) 92 : C WRITE(IOIMP,*) ' 2) compdd' 93 : C WRITE (IOIMP,2019) (compdd(I),I=1,compdd(/2)) 94 : * 95 : * Maintenant, on construit les CORES1 et CORES2 qui nous 96 : * permettent de savoir quels sont les termes diagonaux 97 : * dans les sous-matrices 98 : * 99 : SEGINI CORES1 100 : DO IRIG=1,NRIGEL 101 : DESCR=IRIGEL(3,IRIG) 102 : SEGACT DESCR 103 : * Tout d'abord le segment trav2 qui permet d'avoir 104 : * la correspondance sur les noms d'inconnues 105 : NLIGRP=LISINC(/2) 106 : NLIGRD=LISDUA(/2) 107 : segini trav2 108 : do 410 iligrp=1,nligrp 109 : do 420 incd=1,nbincd 110 : if (compd(incd).eq.lisinc(iligrp)) then 111 : cor1p(iligrp)=incd 112 : goto 410 113 : endif 114 : 420 continue 115 : 410 continue 116 : do 510 iligrd=1,nligrd 117 : do 520 incd=1,nbincd 118 : if (compdd(incd).eq.lisdua(iligrd)) then 119 : cor1d(iligrd)=incd 120 : goto 510 121 : endif 122 : 520 continue 123 : * Si on a pas trouvé dans la liste des duales, on regarde 124 : * dans la liste des primales. 125 : * Ceci est pratique en mécaflu où il arrive qu'on nomme les duales 126 : * pareil que les primales (ex. 'UX' 'UX') 127 : do 525 incd=1,nbincd 128 : if (compd(incd).eq.lisdua(iligrd)) then 129 : cor1d(iligrd)=incd 130 : goto 510 131 : endif 132 : 525 continue 133 : 510 continue 134 : * Voyons si le descripteur est carré, ce qui est un 135 : * cas courant qui simplifie notablement les choses 136 : descar=.false. 137 : if (nligrp.ne.nligrd) goto 600 138 : do iligr=1,nligrp 139 : if (cor1p(iligr).ne.cor1d(iligr)) goto 600 140 : if (noelep(iligr).ne.noeled(iligr)) goto 600 141 : enddo 142 : descar=.true. 143 : 600 continue 144 : * Si oui, tout est simple, sinon il faut rechercher 145 : * une correspondance primale-duale 146 : if (descar) then 147 : ipcor2(irig)=-1 148 : else 149 : nligp=nligrp 150 : segini cores2 151 : nligp=0 152 : do iligrp=1,nligrp 153 : do iligrd=1,nligrd 154 : if (cor1p(iligrp).eq.cor1d(iligrd).and. 155 : $ noelep(iligrp).eq.noeled(iligrd)) then 156 : nligp=nligp+1 157 : cordia(1,nligp)=iligrp 158 : cordia(2,nligp)=iligrd 159 : endif 160 : enddo 161 : enddo 162 : if (nligp.eq.0) then 163 : segsup cores2 164 : ipcor2(irig)=0 165 : else 166 : segadj cores2 167 : ipcor2(irig)=cores2 168 : endif 169 : endif 170 : itrav2(irig)=trav2 171 : C Write(ioimp,*) 'irig=',irig 172 : C Write(ioimp,*) 'descar=',descar 173 : C if (ipcor2.gt.0) then 174 : C write(ioimp,*) 'cordia(1, )' 175 : C WRITE (IOIMP,2020) (cordia(1,I),I=1,cordia(/2)) 176 : C write(ioimp,*) 'cordia(2, )' 177 : C WRITE (IOIMP,2020) (cordia(2,I),I=1,cordia(/2)) 178 : C else 179 : C WRITE(IOIMP,*) 'ipcor2=',ipcor2 180 : C endif 181 : SEGDES DESCR 182 : ENDDO 183 : * 184 : * Maintenant, on peut construire MTRA1, MTRA2 185 : * (repérage des couples noms d'inconnues primales, 186 : * harmoniques) ainsi que CORES3 187 : * et ICPR (repérage global->local des noeuds) 188 : * ce qui permettra dans une deuxième passe de remplir 189 : * le segment TMTRAV qui sera transformé en CHPOINT 190 : * 191 : SEGINI ICPR 192 : NNNOE=0 193 : SEGINI MTRA1,MTRA2 194 : NNIN=0 195 : DO 100 IRIG=1,NRIGEL 196 : DESCR=IRIGEL(3,IRIG) 197 : trav2=itrav2(irig) 198 : MELEME=IRIGEL(1,IRIG) 199 : cores2=ipcor2(irig) 200 : if (cores2.eq.0) goto 99 201 : SEGACT DESCR 202 : SEGACT MELEME 203 : if (cores2.eq.-1) then 204 : nligp=noelep(/1) 205 : else 206 : nligp=cordia(/2) 207 : endif 208 : segini cores3 209 : do 110 iligp=1,nligp 210 : if (cores2.eq.-1) then 211 : iligrp=iligp 212 : else 213 : iligrp=cordia(1,iligp) 214 : endif 215 : do ic=1,icomp(/1) 216 : if (icomp(ic).eq.cor1p(iligrp).and. 217 : $ mhar(ic).eq.irigel(5,irig)) then 218 : comtra(iligp)=ic 219 : goto 110 220 : endif 221 : enddo 222 : icomp(**)=cor1p(iligrp) 223 : mhar(**)=irigel(5,irig) 224 : nnin=nnin+1 225 : comtra(iligp)=nnin 226 : 110 continue 227 : C WRITE(IOIMP,*) ' 3) comtra' 228 : C WRITE (IOIMP,2020) (comtra(I),I=1,comtra(/1)) 229 : ipcor3(irig)=cores3 230 : * Dans le cas cores2=-1 231 : * tous les noeuds du meleme sont référencés 232 : * sinon, on regarde ceux qui le sont 233 : if (cores2.eq.-1) then 234 : nno=num(/1) 235 : else 236 : nno=nligp 237 : endif 238 : do iel=1,num(/2) 239 : do ino=1,nno 240 : if (cores2.eq.-1) then 241 : ilno=ino 242 : else 243 : ilno=noelep(cordia(1,ino)) 244 : endif 245 : ipt=num(ilno,iel) 246 : if (icpr(ipt).eq.0) then 247 : nnnoe=nnnoe+1 248 : icpr(ipt)=nnnoe 249 : endif 250 : enddo 251 : enddo 252 : 99 continue 253 : SEGDES MELEME 254 : SEGDES DESCR 255 : segsup trav2 256 : 100 CONTINUE 257 : C WRITE(IOIMP,*) 'nnin=',nnin 258 : C WRITE(IOIMP,*) 'nnin2=',icomp(/1) 259 : C WRITE(IOIMP,*) 'nnin3=',mhar(/1) 260 : C WRITE(IOIMP,*) ' 1) icomp' 261 : C WRITE (IOIMP,2020) (icomp(I),I=1,icomp(/1)) 262 : C WRITE(IOIMP,*) ' 2) mhar' 263 : C WRITE (IOIMP,2020) (mhar(I),I=1,mhar(/1)) 264 : C WRITE(IOIMP,*) 'nnnoe=',nnnoe 265 : C WRITE(IOIMP,*) ' 4) icpr' 266 : C WRITE (IOIMP,2020) (icpr(I),I=1,icpr(/1)) 267 : * 268 : * Maintenant, on peut remplir MTRAV 269 : * 270 : SEGINI MTRAV 271 : DO ININ=1,NNIN 272 : INCO(ININ)=compd(ICOMP(ININ)) 273 : NHAR(ININ)=MHAR(ININ) 274 : ENDDO 275 : SEGSUP MTRA1,MTRA2 276 : segsup trav1 277 : C xt=0.D0 278 : DO 1100 IRIG=1,NRIGEL 279 : MELEME=IRIGEL(1,IRIG) 280 : DESCR=IRIGEL(3,IRIG) 281 : XMATRI=IRIGEL(4,IRIG) 282 : cores2=ipcor2(irig) 283 : cores3=ipcor3(irig) 284 : if (cores2.eq.0) goto 1099 285 : SEGACT MELEME 286 : SEGACT DESCR 287 : SEGACT XMATRI 288 : nel=num(/2) 289 : if (cores2.eq.-1) then 290 : nligp=noelep(/1) 291 : else 292 : nligp=cordia(/2) 293 : endif 294 : do iel=1,nel 295 : do iligp=1,nligp 296 : if (cores2.eq.-1) then 297 : iligrp=iligp 298 : iligrd=iligrp 299 : else 300 : iligrp=cordia(1,iligp) 301 : iligrd=cordia(2,iligp) 302 : endif 303 : ININ=comtra(iligp) 304 : ino=noelep(iligrp) 305 : ipt=num(ino,iel) 306 : INNOE=icpr(ipt) 307 : C if (ipt.eq.2.and.lisinc(iligrp).eq.'UX ') then 308 : C write(ioimp,*) 'irig=',irig 309 : C write(ioimp,*) 'iel=',iel 310 : C write(ioimp,*) 'iligrp=',iligrp 311 : C write(ioimp,*) 'iligrd=',iligrd 312 : C xt=xt+RE(iligrd,iligrp,iel) 313 : C WRITE(IOIMP,*) 'RE=',RE(iligrd,iligrp,iel) 314 : C WRITE(IOIMP,*) 'xt=',xt 315 : C WRITE(IOIMP,*) 'inin=',inin 316 : C WRITE(IOIMP,*) 'innoe=',innoe 317 : C endif 318 : IGEO(INNOE)=ipt 319 : IBIN(ININ,INNOE)=1 320 : BB(ININ,INNOE)=BB(ININ,INNOE)+ 321 : $ RE(iligrd,iligrp,iel)*COERIG(IRIG) 322 : enddo 323 : enddo 324 : if (cores2.ne.-1) segsup cores2 325 : segsup cores3 326 : SEGDES XMATRI 327 : SEGDES DESCR 328 : SEGDES MELEME 329 : 1099 continue 330 : 1100 CONTINUE 331 : SEGSUP ICPR 332 : SEGSUP CORES1 333 : SEGDES MRIGID 334 : * 335 : * Création du CHPOINT résultat 336 : * 336 : CALL CRECHP(MTRAV,MCHPOI) 337 : $ 338 : SEGSUP MTRAV 339 : RETURN 340 : 2019 FORMAT (20(2X,A4) ) 341 : 2020 FORMAT (20(2X,I4) ) 342 : 2021 FORMAT (20(2X,L4) ) 343 : 2022 FORMAT(10(1X,1PG12.5)) 344 : * 345 : * End of EXDIAG 346 : * 347 : END 348 : 349 : 350 : 351 : 352 : 353 : 354 : 355 : 356 : 357 : 358 : 359 : 360 : 361 : 362 : 363 : 364 : 365 : 366 : 367 : 368 :
© Cast3M 2003 - All rights reserved.
Disclaimer