1 : C RLEXCA SOURCE CHAT 05/01/13 03:02:11 5004 2 : SUBROUTINE RLEXCA(MELLIM,MLELSC,MLESBC,MLRDIS, 3 : & MLESCF,MATCOE) 4 : C************************************************************************ 5 : C 6 : C PROJET : CASTEM 2000 7 : C 8 : C NOM : RLEXCA 9 : C 10 : C DESCRIPTION : Appellé par GRADIA 11 : C 12 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) 13 : C 14 : C AUTEUR : A. BECCANTINI 15 : C 16 : C************************************************************************ 17 : C 18 : C Inputs: 19 : C 20 : C MELLIM : MELEME spg boundary conditions 21 : C 22 : C MLELSC : neighbors of SOMMETS 23 : C 24 : C MLESBC : neighbors of SOMMETS belonging to the boundary 25 : C 26 : C MLRDIS : list of distances of a boundary SOMMET and its neighbors 27 : C 28 : C Output : 29 : C 30 : C MLESCF : final list of SOMMET points and their CENTRE neighbors 31 : C 32 : C MATCOE : coefficient to compute the reconstruction 33 : C 34 : C 35 : C**** Variables de COOPTIO 36 : C 37 : C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR 38 : C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ 39 : C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES 40 : C & ,IECHO, IIMPI, IOSPI 41 : C & ,IDIM 42 : CC & ,MCOORD 43 : C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE 44 : C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU 45 : C & ,NORINC,NORVAL,NORIND,NORVAD 46 : C & ,NUCROU, IPSAUV, IFICLE, IPREFI 47 : CC 48 : IMPLICIT INTEGER(I-N) 49 : -INC CCOPTIO 50 : -INC SMCOORD 51 : -INC SMELEME 52 : POINTEUR MELLIM.MELEME 53 : C 54 : INTEGER NBL, NBTPOI 55 : SEGMENT MLELEM 56 : INTEGER INDEX(NBL+1) 57 : INTEGER LESPOI(NBTPOI) 58 : ENDSEGMENT 59 : POINTEUR MLELSC.MLELEM, MLESBC.MLELEM, MLESCF.MLELEM 60 : C 61 : INTEGER JG 62 : -INC SMLREEL 63 : POINTEUR MLRDIS.MLREEL, MLRSIG.MLREEL, MLRTRA.MLREEL 64 : -INC SMLENTI 65 : POINTEUR MLELIM.MLENTI, MLESOM.MLENTI 66 : C 67 : INTEGER N1,N2 68 : SEGMENT MATRIX 69 : REAL*8 MAT(N1,N2) 70 : ENDSEGMENT 71 : POINTEUR MATA.MATRIX,MATU.MATRIX,MATV.MATRIX 72 : & ,MATCOE.MATRIX 73 : & ,MTAA.MATRIX, MINTAA.MATRIX,MLRIN1.MLREEL,MLRIN2.MLREEL 74 : & ,MLRIN3.MLREEL 75 : C 76 : INTEGER NTP,NBTPO1,NBL1,ISOMM,IPOS,NGS,NLS,IPOS1,NVMAX,IPOS2 77 : & ,IPOS3,NVOI,NVMIN,NSOMM,I1,I2,I3,I4,NGC,IERSVD,NGS1,NVOIS1 78 : & ,NLLIM,IPCOOR,IPVOI,ICEL,NVOI0,IERR0,IVOI 79 : Real*8 XC,YC,ZC,SMAX,SMIN,DIST2,DIST21,ERRTOL,XS,YS,ZS,ERRTO1 80 : PARAMETER(ERRTOL=1.0D-15,ERRTO1=1.0D-7) 81 : LOGICAL LOGSVD 82 : C 83 : NTP=MCOORD.XCOOR(/1)/(IDIM+1) 84 : C 85 : C**** Le MELEME support de CL 86 : C 87 : IF(MELLIM.NE.0) THEN 87 : CALL KRIPAD(MELLIM,MLELIM) 88 : C 89 : C******* En KRIPAD 90 : C SEGACT MELLIM 91 : C SEGINI MLELIM 92 : C 93 : SEGDES MELLIM 94 : ELSE 95 : JG=NTP 96 : SEGINI MLELIM 97 : ENDIF 98 : C 99 : SEGACT MLESBC 100 : NBL1=MLESBC.INDEX(/1)-1 101 : NBTPO1=MLESBC.LESPOI(/1) 102 : JG=NTP 103 : SEGINI MLESOM 104 : DO I1 = 1, NBL1, 1 105 : IPOS=MLESBC.INDEX(I1) 106 : NGS=MLESBC.LESPOI(IPOS) 107 : MLESOM.LECT(NGS)=I1 108 : ENDDO 109 : C 110 : SEGACT MLELSC 111 : NBL=MLELSC.INDEX(/1)-1 112 : NSOMM=NBL 113 : NBTPOI=MLELSC.LESPOI(/1) 114 : NBTPOI=NBTPOI+NBTPO1-NBL1 115 : C 116 : C**** MLESCF 117 : C Structure Sommet-Centres Finale 118 : C NBTPOI(MLESCF) <= NBTPOI (SEGADJ à la fin) 119 : C 120 : SEGINI MLESCF 121 : C 122 : C**** La matrice de coeff. 123 : C 124 : C MATCOE(*,I) = coefficients concernants MLESCF.LESPOI(I) 125 : C 126 : C NVMIN : nombre minimum de voisins que un sommet 127 : C doit avoir pour une reconstruction lineaire 128 : C exacte 129 : C N2(MATCOE) <= N2 (SEGADJ à la fin) 130 : C 131 : NVMIN=IDIM+1 132 : N1=NVMIN 133 : N2=NBTPOI 134 : SEGINI MATCOE 135 : SEGACT MLRDIS 136 : C 137 : C**** Definition des matrices pour calculer MATCOE 138 : C Les dimensions 139 : C 140 : IPOS1=MLELSC.INDEX(1) 141 : NVMAX=0 142 : DO ISOMM = 1, NSOMM, 1 143 : IPOS=IPOS1 144 : IPOS1=MLELSC.INDEX(ISOMM+1) 145 : NVOI=IPOS1-IPOS-1 146 : NGS=MLELSC.LESPOI(IPOS) 147 : NLS=MLESOM.LECT(NGS) 148 : IF(NLS .GT. 0)THEN 149 : IPOS2=MLESBC.INDEX(NLS) 150 : IPOS3=MLESBC.INDEX(NLS+1) 151 : NVOI=NVOI+IPOS3-IPOS2-1 152 : ENDIF 153 : NVMAX=MAX(NVMAX,NVOI) 154 : ENDDO 155 : C 156 : C**** MATA = matrice à inverser (NVMAX,IDIM+1) 157 : C MATU = matice des vectors propres singuliers droite de MATA 158 : C (NVMAX,IDIM+1) 159 : C MATV = matrice des vectors propres singuliers gauche de MATA 160 : C (IDIM+1,IDIM+1) 161 : C Mais en INVSVD a dimension NVMAX,IDIM+1 162 : C MLRSIG =liste des valeurs singuliers de MATA 163 : C 164 : C N.B. MATA = MATU MATSIG t(MATV) 165 : C Si MATA non singulier 166 : C inv(MATA) = MATV inv(MATSIG) t(MATU) 167 : C 168 : N1=NVMAX 169 : N2=NVMIN 170 : SEGINI MATA 171 : SEGINI MATU 172 : SEGINI MATV 173 : JG=NVMIN 174 : SEGINI MLRSIG 175 : SEGINI MLRTRA 176 : C MLRTRA segment de travail de la subroutine invsvd 177 : C 178 : C 179 : C MTAA : [t(MATA).MATA] 180 : C MINTAA : [inve(t(MATA) MATA)] 181 : C MLRIN1,2,3 : "temporary vectors" 182 : C 183 : N1=NVMIN 184 : N2=NVMIN 185 : SEGINI MTAA 186 : SEGINI MINTAA 187 : JG=NVMIN 188 : SEGINI MLRIN1 189 : SEGINI MLRIN2 190 : SEGINI MLRIN3 191 : C 192 : MLESCF.INDEX(1)=1 193 : NBTPOI=0 194 : DO ISOMM = 1, NSOMM, 1 195 : IPOS=MLESCF.INDEX(ISOMM) 196 : IPOS1=MLELSC.INDEX(ISOMM) 197 : NGS=MLELSC.LESPOI(IPOS1) 198 : IPCOOR=(NGS-1)*(IDIM+1)+1 199 : XS=MCOORD.XCOOR(IPCOOR) 200 : YS=MCOORD.XCOOR(IPCOOR+1) 201 : IF(IDIM .EQ. 3) ZS=MCOORD.XCOOR(IPCOOR+2) 202 : MLESCF.LESPOI(IPOS)=NGS 203 : NLLIM=MLELIM.LECT(NGS) 204 : IF(NLLIM.GT.0)THEN 205 : C 206 : C********** NGS est un noeud du MELEME des conditions 207 : C limites 208 : C 209 : MLESCF.INDEX(ISOMM+1)=IPOS+1 210 : MATCOE.MAT(1,IPOS)=1.0D0 211 : DO I2 = 2, NVMIN, 1 212 : MATCOE.MAT(I2,IPOS)=0.0D0 213 : ENDDO 214 : NBTPOI=NBTPOI+1 215 : C 216 : C********** On passe au sommet suivant (i.e. ISOMM->ISOMM+1) 217 : C 218 : ELSE 219 : DO I2 = 1, NVMIN, 1 220 : MATCOE.MAT(I2,IPOS)=0.0D0 221 : ENDDO 222 : NLS=MLESOM.LECT(NGS) 223 : NVOI=MLELSC.INDEX(ISOMM+1)-IPOS1-1 224 : NVOI0=NVOI 225 : C NVOI0 = nombre de centre voisins directe! 226 : DO I2=1,NVOI,1 227 : NGC=MLELSC.LESPOI(IPOS1+I2) 228 : MLESCF.LESPOI(IPOS+I2)=NGC 229 : IPCOOR=(NGC-1)*(IDIM+1)+1 230 : XC=MCOORD.XCOOR(IPCOOR) 231 : YC=MCOORD.XCOOR(IPCOOR+1) 232 : MATA.MAT(I2,1)=1 233 : MATA.MAT(I2,2)=XC-XS 234 : MATA.MAT(I2,3)=YC-YS 235 : IF(IDIM.EQ.3)THEN 236 : ZC=MCOORD.XCOOR(IPCOOR+2) 237 : MATA.MAT(I2,4)=ZC-ZS 238 : ENDIF 239 : ENDDO 240 : IF(NLS.EQ.0)THEN 241 : C 242 : C************* NGS est un noeud À l'interieur du domaine 243 : C 244 : LOGSVD=.TRUE. 244 : CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT, 245 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT,IERSVD, 246 : & MLRTRA.PROG) 247 : IF(IERSVD.NE.0)THEN 248 : C INVSVD did not work 249 : LOGSVD=.FALSE. 250 : ELSE 251 : SMAX=0.0D0 252 : DO I2=1,NVMIN,1 253 : SMAX=MAX(SMAX,MLRSIG.PROG(I2)) 254 : ENDDO 255 : SMIN=SMAX 256 : DO I2=1,NVMIN,1 257 : SMIN=MIN(SMIN,MLRSIG.PROG(I2)) 258 : ENDDO 259 : ENDIF 260 : IF((SMIN/SMAX) .LT. ERRTOL)THEN 261 : C INVSVD did not work 262 : LOGSVD=.FALSE. 263 : ENDIF 264 : ELSE 265 : C 266 : C************* NGS est un noeud sur le bord mais pas un noeud des 267 : C conditions limites 268 : C 269 : IPVOI=MLESBC.INDEX(NLS) 270 : NGS1=MLESBC.LESPOI(IPVOI) 271 : NVOIS1=MLESBC.INDEX(NLS+1)-IPVOI-1 272 : IF((NVOIS1.LE.0).OR.(NGS .NE. NGS1))THEN 273 : WRITE(IOIMP,*) 'Subroutine rlexca.eso.' 273 : CALL ERREUR(5) 274 : GOTO 9999 275 : ENDIF 276 : DIST2=MLRDIS.PROG(IPVOI+1)*1.10D0 277 : ICEL=1 278 : C 279 : C************* Repeat until 280 : C 281 : 10 CONTINUE 282 : NGC=MLESBC.LESPOI(IPVOI+ICEL) 283 : DIST21=MLRDIS.PROG(IPVOI+ICEL) 284 : IF((DIST21.LT.DIST2) .OR. (NVOI .LT. NVMIN))THEN 285 : NVOI=NVOI+1 286 : MLESCF.LESPOI(IPOS+NVOI)=NGC 287 : IPCOOR=(NGC-1)*(IDIM+1)+1 288 : XC=MCOORD.XCOOR(IPCOOR) 289 : YC=MCOORD.XCOOR(IPCOOR+1) 290 : MATA.MAT(NVOI,1)=1 291 : MATA.MAT(NVOI,2)=XC-XS 292 : MATA.MAT(NVOI,3)=YC-YS 293 : IF(IDIM.EQ.3)THEN 294 : ZC=MCOORD.XCOOR(IPCOOR+2) 295 : MATA.MAT(NVOI,4)=ZC-ZS 296 : ENDIF 297 : ICEL=ICEL+1 298 : DIST2=DIST21*1.10D0 299 : IF(ICEL .LE. NVOIS1) GOTO 10 300 : ENDIF 301 : LOGSVD=.TRUE. 301 : CALL INVSVD( NVMAX, NVOI, NVMIN, MATA.MAT, 302 : & MLRSIG.PROG,.TRUE.,MATU.MAT,.TRUE.,MATV.MAT 303 : & ,IERSVD,MLRTRA.PROG) 304 : IF(IERSVD.NE.0)THEN 305 : SMIN=0.0D0 306 : SMAX=1.0D0 307 : ELSE 308 : SMAX=0.0D0 309 : DO I2=1,NVMIN,1 310 : SMAX=MAX(SMAX,MLRSIG.PROG(I2)) 311 : ENDDO 312 : SMIN=SMAX 313 : DO I2=1,NVMIN,1 314 : SMIN=MIN(SMIN,MLRSIG.PROG(I2)) 315 : ENDDO 316 : ENDIF 317 : IF((SMIN/SMAX) .LT. ERRTO1)THEN 318 : C 319 : C**************** We add other neighbors 320 : C 321 : ICEL=NVOI-NVOI0 322 : IF(ICEL.LT.NVOIS1)THEN 323 : ICEL=ICEL+1 324 : NGC=MLESBC.LESPOI(IPVOI+ICEL) 325 : DIST2=MLRDIS.PROG(IPVOI+ICEL)*1.10D0 326 : GOTO 10 327 : ELSE 328 : IF((SMIN/SMAX) .LT. (ERRTOL))THEN 329 : LOGSVD=.FALSE. 330 : ENDIF 331 : ENDIF 332 : ENDIF 333 : C 334 : C************* Fin repeat until 335 : C 336 : ENDIF 337 : C 338 : C********** We have two different cases: 339 : C LOGSVD = .TRUE. -> we have the coeff. we are looking for 340 : C LOGSVD = .FALSE. -> we have not 341 : C 342 : IF(LOGSVD)THEN 343 : DO I4=1,NVMIN,1 344 : DO I3=1,NVOI,1 345 : IPVOI=IPOS+I3 346 : DO I2=1,NVMIN,1 347 : MATCOE.MAT(I2,IPVOI)=MATCOE.MAT(I2,IPVOI)+ 348 : & (MATV.MAT(I2,I4)*MATU.MAT(I3,I4)/ 349 : & MLRSIG.PROG(I4)) 350 : ENDDO 351 : ENDDO 352 : ENDDO 353 : ELSE 354 : WRITE (IOIMP,*) 'rlexca.eso' 355 : C 22 0 356 : C Opération malvenue. Résultat douteux 356 : CALL ERREUR(22) 357 : C 358 : C************* INVSVD does not worked 359 : C For each neighbor k we have to compute the solution 360 : C of 361 : C 362 : C t(MATA) MATA x = t(MATA) * b 363 : C 364 : C where b= \sum_l e_l \delta(k,l) = e_k 365 : C 366 : C To do that, we compute 367 : C 368 : C X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 369 : C 370 : C X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 371 : C [t(MATA) MATA] X_0] 372 : C 373 : C 374 : C************* We compute [t(MATA) MATA] 375 : C We store it in the upper triangle of MTAA(NVMIN,NVMIN) 376 : C 377 : DO I1=1,NVMIN,1 378 : DO I2=I1,NVMIN,1 379 : MTAA.MAT(I1,I2)=0.0D0 380 : ENDDO 381 : ENDDO 382 : C 383 : DO I1=1,NVMIN,1 384 : DO I2=I1,NVMIN,1 385 : DO I3=1,NVOI,1 386 : MTAA.MAT(I1,I2)=MTAA.MAT(I1,I2)+ 387 : & (MATA.MAT(I3,I1)*MATA.MAT(I3,I2)) 388 : ENDDO 389 : ENDDO 390 : ENDDO 391 : C 392 : C************* We compute [inve(t(MATA) MATA)] 393 : C CHOLDC stores it in the upper trianle of MINTAA(NVMIN,NVMIN) 394 : C 394 : CALL CHOLDC(NVMIN,NVMIN,MTAA.MAT,MLRIN1.PROG,MINTAA.MAT, 395 : & MLRIN2.PROG,IERR0) 396 : IF(IERR0.NE.0)THEN 397 : WRITE(IOIMP,*) 'subroutine rlexca.eso.' 398 : C 26 2 399 : C Tache impossible. Probablement données erronées 399 : CALL ERREUR(26) 400 : GOTO 9999 401 : ENDIF 402 : C 403 : C************* We complete MTAA and MINTAA 404 : C 405 : DO I1=1,NVMIN,1 406 : DO I2=I1+1,NVMIN,1 407 : MINTAA.MAT(I2,I1)=MINTAA.MAT(I1,I2) 408 : MTAA.MAT(I2,I1)=MTAA.MAT(I1,I2) 409 : ENDDO 410 : ENDDO 411 : C 412 : DO IVOI=1,NVOI,1 413 : C 414 : C**************** We compute [t(MATA) . b] and we store it in MLRIN1.PROG 415 : C 416 : DO I1=1,NVMIN,1 417 : MLRIN1.PROG(I1)=MATA.MAT(IVOI,I1) 418 : MLRIN2.PROG(I1)=0.0D0 419 : MLRIN3.PROG(I1)=0.0D0 420 : ENDDO 421 : C 422 : C**************** X_0 = [inve(t(MATA) MATA)] [t(MATA) * b] 423 : C X_0(i1) into MLRIN2.PROG(I1) 424 : C 425 : DO I2=1,NVMIN,1 426 : DO I1=1,NVMIN,1 427 : MLRIN2.PROG(I1)=MLRIN2.PROG(I1)+ 428 : & (MINTAA.MAT(I1,I2)*MLRIN1.PROG(I2)) 429 : ENDDO 430 : ENDDO 431 : C 432 : C**************** X_1 = X_0 + [inve(t(MATA) MATA)] [t(MATA) * b - 433 : C [t(MATA) MATA] X_0] 434 : C 435 : C [t(MATA) MATA] X_0 into MLRIN3.PROG 436 : C 437 : DO I2=1,NVMIN,1 438 : DO I1=1,NVMIN,1 439 : MLRIN3.PROG(I1)=MLRIN3.PROG(I1)+ 440 : & (MTAA.MAT(I1,I2)*MLRIN2.PROG(I2)) 441 : ENDDO 442 : ENDDO 443 : C 444 : C**************** Now we have 445 : C [t(MATA) . b] in MLRIN1.PROG 446 : C X_0(i1) in MLRIN2.PROG(I1) 447 : C [t(MATA) MATA] X_0 in MLRIN3.PROG 448 : C 449 : C X_1(i1) in MLRCOE.MAT(i1,IVOI) 450 : C 451 : IPVOI=IPOS+IVOI 452 : DO I1=1,NVMIN,1 453 : DO I2=1,NVMIN,1 454 : MATCOE.MAT(I1,IPVOI)=MATCOE.MAT(I1,IPVOI)+ 455 : & (MINTAA.MAT(I1,I2)* 456 : & (MLRIN1.PROG(I2)-MLRIN3.PROG(I2))) 457 : ENDDO 458 : MATCOE.MAT(I1,IPVOI)=MATCOE.MAT(I1,IPVOI)+ 459 : & MLRIN2.PROG(I1) 460 : ENDDO 461 : ENDDO 462 : ENDIF 463 : CC TEST 464 : C if(isomm .eq. (nsomm - 1))then 465 : C write(*,*) 'ngs =', NGS 466 : C write(*,*) 'invide',LOGSVD 467 : C write(*,*) 'coeff(1) =',(matcoe.mat(1,ipos+ivoi),ivoi=1 468 : C $ ,nvoi,1) 469 : C write(*,*) 'coeff(2) =',(matcoe.mat(2,ipos+ivoi),ivoi=1 470 : C $ ,nvoi,1) 471 : C write(*,*) 'coeff(3) =',(matcoe.mat(3,ipos+ivoi),ivoi=1 472 : C $ ,nvoi,1) 473 : C if(idim.eq.3) write(*,*) 'coeff(3)=' 474 : C & ,(matcoe.mat(4,ipos+ivoi),ivoi=1,nvoi,1) 475 : C dist2=0.0D0 476 : C do ivoi=1,nvoi,1 477 : C dist2=dist2+matcoe.mat(1,ipos+ivoi) 478 : C enddo 479 : C write(*,*) 'sum=',dist2 480 : CC 480 : C call erreur(21) 481 : C goto 9999 482 : C endif 483 : C 484 : C********** If we are here, we have the coeff. in the 'SOMMET' frame. 485 : C 486 : C 487 : NBTPOI=NBTPOI+NVOI+1 488 : MLESCF.INDEX(ISOMM+1)=IPOS+NVOI+1 489 : ENDIF 490 : ENDDO 491 : C 492 : SEGADJ MLESCF 493 : SEGDES MLESCF 494 : C 495 : N2=NBTPOI 496 : SEGADJ MATCOE 497 : SEGDES MATCOE 498 : C 499 : IF(MELLIM .NE. 0) SEGDES MELLIM 500 : SEGSUP MLELIM 501 : SEGSUP MLESBC 502 : SEGSUP MLESOM 503 : SEGSUP MLELSC 504 : SEGSUP MLRDIS 505 : SEGSUP MATA 506 : SEGSUP MATU 507 : SEGSUP MATV 508 : SEGSUP MLRSIG 509 : SEGSUP MLRTRA 510 : C 511 : SEGSUP MTAA 512 : SEGSUP MINTAA 513 : SEGSUP MLRIN1 514 : SEGSUP MLRIN2 515 : SEGSUP MLRIN3 516 : C 517 : 9999 RETURN 518 : END 519 : 520 : 521 : 522 : 523 : 524 : 525 : 526 : 527 : 528 : 529 :
© Cast3M 2003 - All rights reserved.
Disclaimer