OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
matrix.F
Go to the documentation of this file.
1Copyright> OpenRadioss
2Copyright> Copyright (C) 1986-2025 Altair Engineering Inc.
3Copyright>
4Copyright> This program is free software: you can redistribute it and/or modify
5Copyright> it under the terms of the GNU Affero General Public License as published by
6Copyright> the Free Software Foundation, either version 3 of the License, or
7Copyright> (at your option) any later version.
8Copyright>
9Copyright> This program is distributed in the hope that it will be useful,
10Copyright> but WITHOUT ANY WARRANTY; without even the implied warranty of
11Copyright> MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12Copyright> GNU Affero General Public License for more details.
13Copyright>
14Copyright> You should have received a copy of the GNU Affero General Public License
15Copyright> along with this program. If not, see <https://www.gnu.org/licenses/>.
16Copyright>
17Copyright>
18Copyright> Commercial Alternative: Altair Radioss Software
19Copyright>
20Copyright> As an alternative to this open-source version, Altair also offers Altair Radioss
21Copyright> software under a commercial license. Contact Altair to discuss further if the
22Copyright> commercial version may interest you: https://www.altair.com/radioss/.
23!||====================================================================
24!|| jacobiew ../starter/source/materials/tools/matrix.F
25!||--- called by ------------------------------------------------------
26!|| sigeps38 ../starter/source/materials/mat/mat038/sigeps38.F
27!||====================================================================
28 SUBROUTINE jacobiew(A,N,EW,EV,NROT)
29C-------------------------------------------------------KK141189-
30C COMPUTATION OF ALL EIGENVALUES AND EIGENVECTORS OF A SYMMETRIC
31C MATRIX A BY THE JACOBI ALGORITHM
32C
33C A(N,N) EIGENWERTPROBLEM
34C N DIMENSION OF A
35C EW(N) EIGENVALUES
36C EV(N,N) EIGENVEKTORS
37C NROT NUMBER OF ROTATIONS
38C MAXA MAXIMUM ELEMENT OF A
39C-----------------------------------------------
40C I m p l i c i t T y p e s
41C-----------------------------------------------
42#include "implicit_f.inc"
43 INTEGER NN
44 parameter(nn=9)
45
46 INTEGER N,NROT
47 my_real
48 . a(n,n), ew(n), ev(n,n)
49 . , b(nn), z(nn)
50 INTEGER IZ,IS,ITER,J
51 my_real
52 . sumrs,eps,g,h,t,c,s,tau,theta
53C----------------------------------------------------------------
54 DO 130 iz=1,n
55 DO 120 is=1,n
56C PRACTICAL FOR RADIOSS (PASS ONLY LOWER DIAGONAL MATRIX)
57 IF(iz>is) a(is,iz) = a(iz,is)
58 ev(iz,is)=zero
59 120 CONTINUE
60 b(iz)=a(iz,iz)
61 ew(iz)=b(iz)
62 z(iz)=0.
63C ... EV(IZ,IZ)=1.0
64 ev(iz,iz)=zero
65 130 CONTINUE
66
67 nrot=0
68
69C START ITERATION
70
71 DO 240 iter = 1,50
72 sumrs = zero
73
74C SUM OF THE OFF DIAGONALS
75 DO 150 iz=1,n-1
76 DO 140 is=iz+1,n
77 sumrs=sumrs+abs(a(iz,is))
78 140 CONTINUE
79 150 CONTINUE
80
81 IF (sumrs ==zero) GOTO 9000
82 IF (iter > 4) THEN
83 eps = zero
84 ELSE
85 eps = one_fifth*sumrs/n**2
86 ENDIF
87
88 DO 220 iz=1,n-1
89 DO 210 is=iz+1,n
90 g = 100. * abs(a(iz,is))
91 IF (iter>4 .AND. abs(ew(iz))+g==abs(ew(iz))
92 & .AND. abs(ew(is))+g==abs(ew(is))) THEN
93 a(iz,is)=zero
94 ELSE IF (abs(a(iz,is)) > eps) THEN
95 h = ew(is)-ew(iz)
96 IF (abs(h)+g==abs(h)) THEN
97 t = a(iz,is)/h
98 ELSE
99 theta = half*h/a(iz,is)
100 t=one/(abs(theta)+sqrt(one+theta**2))
101 IF (theta < zero) t=-t
102 ENDIF
103 c=one/sqrt(one+t**2)
104 s=t*c
105 tau=s/(one+c)
106 h=t*a(iz,is)
107 z(iz)=z(iz)-h
108 z(is)=z(is)+h
109 ew(iz)=ew(iz)-h
110 ew(is)=ew(is)+h
111 a(iz,is)=zero
112 DO 160 j=1,iz-1
113 g=a(j,iz)
114 h=a(j,is)
115 a(j,iz)=g-s*(h+g*tau)
116 a(j,is)=h+s*(g-h*tau)
117 160 CONTINUE
118 DO 170 j=iz+1,is-1
119 g=a(iz,j)
120 h=a(j,is)
121 a(iz,j)=g-s*(h+g*tau)
122 a(j,is)=h+s*(g-h*tau)
123 170 CONTINUE
124 DO 180 j=is+1,n
125 g=a(iz,j)
126 h=a(is,j)
127 a(iz,j)=g-s*(h+g*tau)
128 a(is,j)=h+s*(g-h*tau)
129 180 CONTINUE
130 DO 190 j=1,n
131 g=ev(j,iz)
132 h=ev(j,is)
133 ev(j,iz)=g-s*(h+g*tau)
134 ev(j,is)=h+s*(g-h*tau)
135 190 CONTINUE
136 nrot=nrot+1
137 ENDIF
138 210 CONTINUE
139 220 CONTINUE
140
141 DO 230 iz=1,n
142 b(iz)=b(iz)+z(iz)
143 ew(iz)=b(iz)
144 z(iz)=zero
145 230 CONTINUE
146 240 CONTINUE
147
148 9000 CONTINUE
149
150 RETURN
151
152 END
153
154
155!||====================================================================
156!|| dreh ../starter/source/materials/tools/matrix.F
157!||--- called by ------------------------------------------------------
158!|| sigeps38 ../starter/source/materials/mat/mat038/sigeps38.F
159!||====================================================================
160 SUBROUTINE dreh(B,TR,II,JJ,KEN)
161C-------------------------------------------------------KK010587-
162C KEN=0 -> ROTATION OF MATRIX B BY TR = TR(T)*B*TR
163C KEN=1 -> ROTATION OF MATRIX B BY TR(T) = TR*B*TR(T)
164C II AND JJ POINT TO THE SUBMATRIX TO BE ROTATED; DEFAULT II=JJ=1
165C----------------------------------------------------------------
166C-----------------------------------------------
167C I m p l i c i t T y p e s
168C-----------------------------------------------
169#include "implicit_f.inc"
170 my_real
171 . b(3,3),tr(3,3),lk(3,3),x
172 INTEGER KEN,II,JJ
173 INTEGER I,J,K,I1,J1
174C----------------------------------------------------------------
175 IF(ii<=0) ii=1
176 IF(jj<=0) jj=1
177
178 DO 20 i=1,3
179 DO 20 j=1,3
180 j1=jj+j-1
181 x=0.0
182 DO 10 k=1,3
183C IF(TR(K,I)==ZERO) GOTO 10
184 i1=k+ii-1
185C IF(B(I1,J1)==ZERO) GOTO 10
186 IF(ken/=1) x=x+tr(k,i)*b(i1,j1)
187 IF(ken==1) x=x+tr(i,k)*b(i1,j1)
188 10 CONTINUE
189 20 lk(i,j)=x
190
191 DO 40 i=1,3
192 DO 40 j=1,3
193 x=zero
194 DO 30 k=1,3
195C IF(TR(K,J)==ZERO) GOTO 30
196C IF(LK(I,K)==ZERO) GOTO 30
197 IF(ken/=1) x=x+lk(i,k)*tr(k,j)
198 IF(ken==1) x=x+lk(i,k)*tr(j,k)
199 30 CONTINUE
200 40 b(ii+i-1,jj+j-1)=x
201
202 RETURN
203
204 END
205C
206!||====================================================================
207!|| valpvec ../starter/source/materials/tools/matrix.F
208!||--- called by ------------------------------------------------------
209!|| sigeps38 ../starter/source/materials/mat/mat038/sigeps38.F
210!|| sigeps42 ../starter/source/materials/mat/mat042/sigeps42.F
211!|| sigeps90 ../starter/source/materials/mat/mat090/sigeps90.F
212!||--- calls -----------------------------------------------------
213!||====================================================================
214 SUBROUTINE valpvec(SIG,VAL,VEC,NEL)
215C-----------------------------------------------
216C I m p l i c i t T y p e s
217C-----------------------------------------------
218#include "implicit_f.inc"
219C-----------------------------------------------
220C G l o b a l P a r a m e t e r s
221C-----------------------------------------------
222#include "mvsiz_p.inc"
223C-----------------------------------------------
224C D u m m y A r g u m e n t s
225C-----------------------------------------------
226 my_real
227 . sig(6,*), val(3,*), vec(9,*)
228 INTEGER NEL
229C-----------------------------------------------
230C L o c a l V a r i a b l e s
231C-----------------------------------------------
232 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
233 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
234 my_real
235 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
236 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
237 . cc, angp, dd, ftpi, ttpi, strmax,
238 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
239 . vmag, s11,
240 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
241 . a22, a23, a31, a32, a33,
242 . mdemi, xmaxinv, flm
243 REAL FLMIN
244C-----------------------------------------------
245C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
246C FTPI=(4/3)*PI, TTPI=(2/3)*PI
247C
248C DEVIATEUR PRINCIPAL DE CONTRAINTE
249C . . . . . . . . . . . . . . . . . . .
250 mdemi = -half
251 ttpi = acos(mdemi)
252 ftpi = two*ttpi
253C precision minimum dependant du type REAL ou DOUBLE
254 CALL floatmin(cs(1),cs(2),flmin)
255 flm = two*sqrt(flmin)
256 nindex3=0
257 DO nn = 1, nel
258 cs(1) = sig(1,nn)
259 cs(2) = sig(2,nn)
260 cs(3) = sig(3,nn)
261 cs(4) = sig(4,nn)
262 cs(5) = sig(5,nn)
263 cs(6) = sig(6,nn)
264 pr = -(cs(1)+cs(2)+cs(3))* third
265 cs(1) = cs(1) + pr
266 cs(2) = cs(2) + pr
267 cs(3) = cs(3) + pr
268 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
269 & - cs(2)*cs(3) - cs(1)*cs(3)
270 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
271 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
272 norminf(nn) = em10*norminf(nn)
273C cas racine triple
274 aa = max(aaa(nn),norminf(nn),em20)
275C
276 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
277 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
278 & - two*cs(4)*cs(5)*cs(6)
279C
280 cc=-sqrt(twenty7/aa)*bb*half/aa
281 cc= min(cc,one)
282 cc= max(cc,-one)
283 angp=acos(cc) * third
284 dd=two*sqrt(aa * third)
285 str(1,nn)=dd*cos(angp)
286 str(2,nn)=dd*cos(angp+ftpi)
287 str(3,nn)=dd*cos(angp+ttpi)
288C
289 val(1,nn) = str(1,nn)-pr
290 val(2,nn) = str(2,nn)-pr
291 val(3,nn) = str(3,nn)-pr
292C renforcement de precision en compression
293 IF(abs(str(3,nn))>abs(str(1,nn))
294 & .AND.aaa(nn)>norminf(nn)) THEN
295 aa=str(1,nn)
296 str(1,nn)=str(3,nn)
297 str(3,nn)=aa
298 nindex3 = nindex3+1
299 index3(nindex3) = nn
300 ENDIF
301C . . . . . . . . . . .
302C VECTEURS PROPRES
303C . . . . . . . . . . .
304 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
305 tol1(nn)= max(em20,flm*strmax**2)
306 tol2(nn)=flm*strmax/3
307 a(1,1,nn)=cs(1)-str(1,nn)
308 a(2,2,nn)=cs(2)-str(1,nn)
309 a(3,3,nn)=cs(3)-str(1,nn)
310 a(1,2,nn)=cs(4)
311 a(2,1,nn)=cs(4)
312 a(2,3,nn)=cs(5)
313 a(3,2,nn)=cs(5)
314 a(1,3,nn)=cs(6)
315 a(3,1,nn)=cs(6)
316C
317 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
318 . *a(2,2,nn)
319 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
320 . *a(2,3,nn)
321 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
322 . *a(2,1,nn)
323 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
324 . *a(3,2,nn)
325 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
326 . *a(3,3,nn)
327 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
328 . *a(3,1,nn)
329 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
330 . *a(1,2,nn)
331 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
332 . *a(1,3,nn)
333 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
334 . *a(1,1,nn)
335 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
336 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
337 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
338
339 ENDDO
340C
341 nindex1 = 0
342 nindex2 = 0
343 DO nn = 1, nel
344 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
345 IF(xmag(1,nn)==xmaxv(nn)) THEN
346 lmaxv(nn) = 1
347 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
348 lmaxv(nn) = 2
349 ELSE
350 lmaxv(nn) = 3
351 ENDIF
352 IF(aaa(nn)<norminf(nn)) THEN
353 val(1,nn) = sig(1,nn)
354 val(2,nn) = sig(2,nn)
355 val(3,nn) = sig(3,nn)
356 v(1,1,nn) = one
357 v(2,1,nn) = zero
358 v(3,1,nn) = zero
359 v(1,2,nn) = zero
360 v(2,2,nn) = one
361 v(3,2,nn) = zero
362
363 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
364 nindex1 = nindex1 + 1
365 index1(nindex1) = nn
366 ELSE
367 nindex2 = nindex2 + 1
368 index2(nindex2) = nn
369 ENDIF
370 ENDDO
371C
372C #include "vectorize.inc"
373 DO n = 1, nindex1
374 nn = index1(n)
375 lmax = lmaxv(nn)
376 xmaxinv = one/xmaxv(nn)
377 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
378 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
379 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
380 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
381 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
382 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
383C
384 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
385 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
386 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
387 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
388 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
389 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
390 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
391 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
392 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
393 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
394 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
395 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
396C
397 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
398 ENDDO
399C
400C #include "vectorize.inc"
401 DO n = 1, nindex1
402 nn = index1(n)
403 IF(xmag(1,nn)==xmaxv(nn)) THEN
404 lmaxv(nn) = 1
405 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
406 lmaxv(nn) = 2
407 ELSE
408 lmaxv(nn) = 3
409 ENDIF
410C
411 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
412 lmax = lmaxv(nn)
413 IF(xmaxv(nn)>tol2(nn))THEN
414 xmaxinv = one/xmaxv(nn)
415 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
416 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
417 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
418 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
419 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
420 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
421 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
422 v(1,2,nn)=v(1,2,nn)*vmag
423 v(2,2,nn)=v(2,2,nn)*vmag
424 v(3,2,nn)=v(3,2,nn)*vmag
425 ELSEIF(vmag>tol2(nn))THEN
426 v(1,2,nn)=-v(2,1,nn)/vmag
427 v(2,2,nn)=v(1,1,nn)/vmag
428 v(3,2,nn)=zero
429 ELSE
430 v(1,2,nn)=one
431 v(2,2,nn)=zero
432 v(3,2,nn)=zero
433 ENDIF
434 ENDDO
435C . . . . . . . . . . . . .
436C SOLUTION DOUBLE
437C . . . . . . . . . . . . .
438 DO n = 1, nindex2
439 nn = index2(n)
440 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
441 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
442 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
443C
444 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
445 ENDDO
446C
447C #include "vectorize.inc"
448 DO n = 1, nindex2
449 nn = index2(n)
450 IF(xmag(1,nn)==xmaxv(nn)) THEN
451 lmaxv(nn) = 1
452 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
453 lmaxv(nn) = 2
454 ELSE
455 lmaxv(nn) = 3
456 ENDIF
457C
458 lmax = lmaxv(nn)
459 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
460 & <tol2(nn))THEN
461 xmaxinv = one/xmaxv(nn)
462 v(1,1,nn)= zero
463 v(2,1,nn)= zero
464 v(3,1,nn)= one
465 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
466 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
467 v(3,2,nn)= zero
468C
469 ELSEIF(xmaxv(nn)>tol2(nn))THEN
470 xmaxinv = one/xmaxv(nn)
471 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
472 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
473 v(3,1,nn)= zero
474 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
475 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
476 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
477 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
478 v(1,2,nn)=v(1,2,nn)*vmag
479 v(2,2,nn)=v(2,2,nn)*vmag
480 v(3,2,nn)=v(3,2,nn)*vmag
481 ELSE
482 v(1,1,nn) = one
483 v(2,1,nn) = zero
484 v(3,1,nn) = zero
485 v(1,2,nn) = zero
486 v(2,2,nn) = one
487 v(3,2,nn) = zero
488 ENDIF
489 ENDDO
490C
491 DO nn = 1, nel
492 vec(1,nn)=v(1,1,nn)
493 vec(2,nn)=v(2,1,nn)
494 vec(3,nn)=v(3,1,nn)
495 vec(4,nn)=v(1,2,nn)
496 vec(5,nn)=v(2,2,nn)
497 vec(6,nn)=v(3,2,nn)
498 vec(7,nn)=vec(2,nn)*vec(6,nn)-vec(3,nn)*vec(5,nn)
499 vec(8,nn)=vec(3,nn)*vec(4,nn)-vec(1,nn)*vec(6,nn)
500 vec(9,nn)=vec(1,nn)*vec(5,nn)-vec(2,nn)*vec(4,nn)
501 ENDDO
502C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
503 DO n = 1, nindex3
504 nn = index3(n)
505C str utilise com tableau temporaire au lieu de scalaires temporaires
506 str(1,nn)=vec(7,nn)
507 str(2,nn)=vec(8,nn)
508 str(3,nn)=vec(9,nn)
509 ENDDO
510 DO n = 1, nindex3
511 nn = index3(n)
512 vec(7,nn)=vec(1,nn)
513 vec(8,nn)=vec(2,nn)
514 vec(9,nn)=vec(3,nn)
515 vec(1,nn)=-str(1,nn)
516 vec(2,nn)=-str(2,nn)
517 vec(3,nn)=-str(3,nn)
518 ENDDO
519C
520 RETURN
521 END
522
523!||====================================================================
524!|| valpvecdp ../starter/source/materials/tools/matrix.F
525!||--- called by ------------------------------------------------------
526!|| sigeps01 ../starter/source/materials/mat/mat001/sigeps01.f90
527!|| sigeps38 ../starter/source/materials/mat/mat038/sigeps38.F
528!|| sigeps42 ../starter/source/materials/mat/mat042/sigeps42.F
529!|| sigeps90 ../starter/source/materials/mat/mat090/sigeps90.F
530!||--- calls -----------------------------------------------------
531!||====================================================================
532 SUBROUTINE valpvecdp(SIG,VAL,VEC,NEL)
533C-----------------------------------------------
534C I m p l i c i t T y p e s
535C-----------------------------------------------
536#include "implicit_f.inc"
537C-----------------------------------------------
538C G l o b a l P a r a m e t e r s
539C-----------------------------------------------
540#include "mvsiz_p.inc"
541C-----------------------------------------------
542C D u m m y A r g u m e n t s
543C-----------------------------------------------
544 my_real
545 . sig(6,*), val(3,*), vec(9,*)
546 INTEGER NEL
547C-----------------------------------------------
548C L o c a l V a r i a b l e s
549C-----------------------------------------------
550 INTEGER I, L, N, NN, LMAX, LMAXV(MVSIZ),NINDEX3, INDEX3(MVSIZ),
551 . NINDEX1, NINDEX2, INDEX1(MVSIZ), INDEX2(MVSIZ)
552 double precision
553 . cs(6), str(3,mvsiz), a(3,3,mvsiz), v(3,3,mvsiz), b(3,3,mvsiz),
554 . xmag(3,mvsiz), pr, aa, bb, aaa(mvsiz),
555 . cc, angp, dd, ftpi, ttpi, strmax,
556 . tol1(mvsiz), tol2(mvsiz), xmaxv(mvsiz), norminf(mvsiz),
557 . vmag, s11,
558 . s21, s31, s12, s22, s32, s13, s23, s33, a11, a12, a13, a21,
559 . a22, a23, a31, a32, a33,
560 . mdemi, xmaxinv, flm,
561 . valdp(3,mvsiz),vecdp(9,mvsiz)
562 REAL FLMIN
563C-----------------------------------------------
564C DATA FTPI,TTPI / 4.188790205, 2.094395102 /
565C FTPI=(4/3)*PI, TTPI=(2/3)*PI
566C
567C DEVIATEUR PRINCIPAL DE CONTRAINTE
568C . . . . . . . . . . . . . . . . . . .
569 mdemi = -half
570 ttpi = acos(mdemi)
571 ftpi = two*ttpi
572C precision minimum dependant du type REAL ou DOUBLE
573 CALL floatmin(cs(1),cs(2),flmin)
574 flm = two*sqrt(flmin)
575 nindex3=0
576 DO nn = 1, nel
577 cs(1) = sig(1,nn)
578 cs(2) = sig(2,nn)
579 cs(3) = sig(3,nn)
580 cs(4) = sig(4,nn)
581 cs(5) = sig(5,nn)
582 cs(6) = sig(6,nn)
583 pr = -(cs(1)+cs(2)+cs(3)) * third
584 cs(1) = cs(1) + pr
585 cs(2) = cs(2) + pr
586 cs(3) = cs(3) + pr
587 aaa(nn)=cs(4)**2 + cs(5)**2 + cs(6)**2 - cs(1)*cs(2)
588 & - cs(2)*cs(3) - cs(1)*cs(3)
589 norminf(nn) = max(abs(cs(1)),abs(cs(2)),abs(cs(3)),
590 & abs(cs(4)),abs(cs(5)),abs(cs(6)))
591 norminf(nn) = em10*norminf(nn)
592C cas racine triple
593 aa = max(aaa(nn),norminf(nn),em20)
594C
595 bb=cs(1)*cs(5)**2 + cs(2)*cs(6)**2
596 & + cs(3)*cs(4)**2 - cs(1)*cs(2)*cs(3)
597 & - two*cs(4)*cs(5)*cs(6)
598C
599 cc=-sqrt(twenty7/aa)*bb*half/aa
600 cc= min(cc,one)
601 cc= max(cc,-one)
602 angp=acos(cc) * third
603 dd=two*sqrt(aa * third)
604 str(1,nn)=dd*cos(angp)
605 str(2,nn)=dd*cos(angp+ftpi)
606 str(3,nn)=dd*cos(angp+ttpi)
607C
608 valdp(1,nn) = str(1,nn)-pr
609 valdp(2,nn) = str(2,nn)-pr
610 valdp(3,nn) = str(3,nn)-pr
611C renforcement de precision en compression simple----
612 IF(abs(str(3,nn))>abs(str(1,nn))
613 & .AND.aaa(nn)>norminf(nn)) THEN
614 aa=str(1,nn)
615 str(1,nn)=str(3,nn)
616 str(3,nn)=aa
617 nindex3 = nindex3+1
618 index3(nindex3) = nn
619 ENDIF
620C . . . . . . . . . . .
621C VECTEURS PROPRES
622C . . . . . . . . . . .
623 strmax= max(abs(str(1,nn)),abs(str(3,nn)))
624 tol1(nn)= max(em20,flm*strmax**2)
625 tol2(nn)=flm*strmax * third
626 a(1,1,nn)=cs(1)-str(1,nn)
627 a(2,2,nn)=cs(2)-str(1,nn)
628 a(3,3,nn)=cs(3)-str(1,nn)
629 a(1,2,nn)=cs(4)
630 a(2,1,nn)=cs(4)
631 a(2,3,nn)=cs(5)
632 a(3,2,nn)=cs(5)
633 a(1,3,nn)=cs(6)
634 a(3,1,nn)=cs(6)
635C
636 b(1,1,nn)=a(2,1,nn)*a(3,2,nn)-a(3,1,nn)
637 . *a(2,2,nn)
638 b(1,2,nn)=a(2,2,nn)*a(3,3,nn)-a(3,2,nn)
639 . *a(2,3,nn)
640 b(1,3,nn)=a(2,3,nn)*a(3,1,nn)-a(3,3,nn)
641 . *a(2,1,nn)
642 b(2,1,nn)=a(3,1,nn)*a(1,2,nn)-a(1,1,nn)
643 . *a(3,2,nn)
644 b(2,2,nn)=a(3,2,nn)*a(1,3,nn)-a(1,2,nn)
645 . *a(3,3,nn)
646 b(2,3,nn)=a(3,3,nn)*a(1,1,nn)-a(1,3,nn)
647 . *a(3,1,nn)
648 b(3,1,nn)=a(1,1,nn)*a(2,2,nn)-a(2,1,nn)
649 . *a(1,2,nn)
650 b(3,2,nn)=a(1,2,nn)*a(2,3,nn)-a(2,2,nn)
651 . *a(1,3,nn)
652 b(3,3,nn)=a(1,3,nn)*a(2,1,nn)-a(2,3,nn)
653 . *a(1,1,nn)
654 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
655 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
656 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
657
658 ENDDO
659C
660 nindex1 = 0
661 nindex2 = 0
662 DO nn = 1, nel
663 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
664 IF(xmag(1,nn)==xmaxv(nn)) THEN
665 lmaxv(nn) = 1
666 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
667 lmaxv(nn) = 2
668 ELSE
669 lmaxv(nn) = 3
670 ENDIF
671 IF(aaa(nn)<norminf(nn)) THEN
672 valdp(1,nn) = sig(1,nn)
673 valdp(2,nn) = sig(2,nn)
674 valdp(3,nn) = sig(3,nn)
675 v(1,1,nn) = one
676 v(2,1,nn) = zero
677 v(3,1,nn) = zero
678 v(1,2,nn) = zero
679 v(2,2,nn) = one
680 v(3,2,nn) = zero
681 ELSEIF(xmaxv(nn)>tol1(nn)) THEN
682 nindex1 = nindex1 + 1
683 index1(nindex1) = nn
684 ELSE
685 nindex2 = nindex2 + 1
686 index2(nindex2) = nn
687 ENDIF
688 ENDDO
689C
690C #include "vectorize.inc"
691 DO n = 1, nindex1
692 nn = index1(n)
693 lmax = lmaxv(nn)
694 xmaxinv = one/xmaxv(nn)
695 v(1,1,nn)=b(1,lmax,nn)*xmaxinv
696 v(2,1,nn)=b(2,lmax,nn)*xmaxinv
697 v(3,1,nn)=b(3,lmax,nn)*xmaxinv
698 a(1,1,nn)=a(1,1,nn)+str(1,nn)-str(3,nn)
699 a(2,2,nn)=a(2,2,nn)+str(1,nn)-str(3,nn)
700 a(3,3,nn)=a(3,3,nn)+str(1,nn)-str(3,nn)
701C
702 b(1,1,nn)=a(2,1,nn)*v(3,1,nn)-a(3,1,nn)*v(2,1,nn)
703 b(1,2,nn)=a(2,2,nn)*v(3,1,nn)-a(3,2,nn)*v(2,1,nn)
704 b(1,3,nn)=a(2,3,nn)*v(3,1,nn)-a(3,3,nn)*v(2,1,nn)
705 b(2,1,nn)=a(3,1,nn)*v(1,1,nn)-a(1,1,nn)*v(3,1,nn)
706 b(2,2,nn)=a(3,2,nn)*v(1,1,nn)-a(1,2,nn)*v(3,1,nn)
707 b(2,3,nn)=a(3,3,nn)*v(1,1,nn)-a(1,3,nn)*v(3,1,nn)
708 b(3,1,nn)=a(1,1,nn)*v(2,1,nn)-a(2,1,nn)*v(1,1,nn)
709 b(3,2,nn)=a(1,2,nn)*v(2,1,nn)-a(2,2,nn)*v(1,1,nn)
710 b(3,3,nn)=a(1,3,nn)*v(2,1,nn)-a(2,3,nn)*v(1,1,nn)
711 xmag(1,nn)=sqrt(b(1,1,nn)**2+b(2,1,nn)**2+b(3,1,nn)**2)
712 xmag(2,nn)=sqrt(b(1,2,nn)**2+b(2,2,nn)**2+b(3,2,nn)**2)
713 xmag(3,nn)=sqrt(b(1,3,nn)**2+b(2,3,nn)**2+b(3,3,nn)**2)
714C
715 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
716 ENDDO
717C
718C #include "vectorize.inc"
719 DO n = 1, nindex1
720 nn = index1(n)
721 IF(xmag(1,nn)==xmaxv(nn)) THEN
722 lmaxv(nn) = 1
723 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
724 lmaxv(nn) = 2
725 ELSE
726 lmaxv(nn) = 3
727 ENDIF
728C
729 vmag=sqrt(v(1,1,nn)**2+v(2,1,nn)**2)
730 lmax = lmaxv(nn)
731 IF(xmaxv(nn)>tol2(nn))THEN
732 xmaxinv = one/xmaxv(nn)
733 v(1,3,nn)=b(1,lmax,nn)*xmaxinv
734 v(2,3,nn)=b(2,lmax,nn)*xmaxinv
735 v(3,3,nn)=b(3,lmax,nn)*xmaxinv
736 v(1,2,nn)=v(2,3,nn)*v(3,1,nn)-v(2,1,nn)*v(3,3,nn)
737 v(2,2,nn)=v(3,3,nn)*v(1,1,nn)-v(3,1,nn)*v(1,3,nn)
738 v(3,2,nn)=v(1,3,nn)*v(2,1,nn)-v(1,1,nn)*v(2,3,nn)
739 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
740 v(1,2,nn)=v(1,2,nn)*vmag
741 v(2,2,nn)=v(2,2,nn)*vmag
742 v(3,2,nn)=v(3,2,nn)*vmag
743 ELSEIF(vmag>tol2(nn))THEN
744 v(1,2,nn)=-v(2,1,nn)/vmag
745 v(2,2,nn)=v(1,1,nn)/vmag
746 v(3,2,nn)=zero
747 ELSE
748 v(1,2,nn)=one
749 v(2,2,nn)=zero
750 v(3,2,nn)=zero
751 ENDIF
752 ENDDO
753C . . . . . . . . . . . . .
754C SOLUTION DOUBLE
755C . . . . . . . . . . . . .
756 DO n = 1, nindex2
757 nn = index2(n)
758 xmag(1,nn)=sqrt(a(1,1,nn)**2+a(2,1,nn)**2)
759 xmag(2,nn)=sqrt(a(1,2,nn)**2+a(2,2,nn)**2)
760 xmag(3,nn)=sqrt(a(1,3,nn)**2+a(2,3,nn)**2)
761C
762 xmaxv(nn)=max(xmag(1,nn),xmag(2,nn),xmag(3,nn))
763 ENDDO
764C
765C #include "vectorize.inc"
766 DO n = 1, nindex2
767 nn = index2(n)
768 IF(xmag(1,nn)==xmaxv(nn)) THEN
769 lmaxv(nn) = 1
770 ELSEIF(xmag(2,nn)==xmaxv(nn)) THEN
771 lmaxv(nn) = 2
772 ELSE
773 lmaxv(nn) = 3
774 ENDIF
775C
776 lmax = lmaxv(nn)
777 IF(max(abs(a(3,1,nn)),abs(a(3,2,nn)),abs(a(3,3,nn)))
778 & <tol2(nn))THEN
779 xmaxinv = one/xmaxv(nn)
780 v(1,1,nn)= zero
781 v(2,1,nn)= zero
782 v(3,1,nn)= one
783 v(1,2,nn)=-a(2,lmax,nn)*xmaxinv
784 v(2,2,nn)= a(1,lmax,nn)*xmaxinv
785 v(3,2,nn)= zero
786C
787 ELSEIF(xmaxv(nn)>tol2(nn))THEN
788 xmaxinv = one/xmaxv(nn)
789 v(1,1,nn)=-a(2,lmax,nn)*xmaxinv
790 v(2,1,nn)= a(1,lmax,nn)*xmaxinv
791 v(3,1,nn)= zero
792 v(1,2,nn)=-a(3,lmax,nn)*v(2,1,nn)
793 v(2,2,nn)= a(3,lmax,nn)*v(1,1,nn)
794 v(3,2,nn)= a(1,lmax,nn)*v(2,1,nn)-a(2,lmax,nn)*v(1,1,nn)
795 vmag=one/sqrt(v(1,2,nn)**2+v(2,2,nn)**2+v(3,2,nn)**2)
796 v(1,2,nn)=v(1,2,nn)*vmag
797 v(2,2,nn)=v(2,2,nn)*vmag
798 v(3,2,nn)=v(3,2,nn)*vmag
799 ELSE
800 v(1,1,nn) = one
801 v(2,1,nn) = zero
802 v(3,1,nn) = zero
803 v(1,2,nn) = zero
804 v(2,2,nn) = one
805 v(3,2,nn) = zero
806 ENDIF
807 ENDDO
808C
809 DO nn = 1, nel
810 vecdp(1,nn)=v(1,1,nn)
811 vecdp(2,nn)=v(2,1,nn)
812 vecdp(3,nn)=v(3,1,nn)
813 vecdp(4,nn)=v(1,2,nn)
814 vecdp(5,nn)=v(2,2,nn)
815 vecdp(6,nn)=v(3,2,nn)
816 vecdp(7,nn)=vecdp(2,nn)*vecdp(6,nn)-vecdp(3,nn)*vecdp(5,nn)
817 vecdp(8,nn)=vecdp(3,nn)*vecdp(4,nn)-vecdp(1,nn)*vecdp(6,nn)
818 vecdp(9,nn)=vecdp(1,nn)*vecdp(5,nn)-vecdp(2,nn)*vecdp(4,nn)
819 ENDDO
820C
821 DO nn = 1, nel
822 val(1,nn)=valdp(1,nn)
823 val(2,nn)=valdp(2,nn)
824 val(3,nn)=valdp(3,nn)
825 vec(1,nn)=vecdp(1,nn)
826 vec(2,nn)=vecdp(2,nn)
827 vec(3,nn)=vecdp(3,nn)
828 vec(4,nn)=vecdp(4,nn)
829 vec(5,nn)=vecdp(5,nn)
830 vec(6,nn)=vecdp(6,nn)
831 vec(7,nn)=vecdp(7,nn)
832 vec(8,nn)=vecdp(8,nn)
833 vec(9,nn)=vecdp(9,nn)
834 ENDDO
835C reecriture pour contourner probleme sur itanium2 comp 9. + latency=16
836 DO n = 1, nindex3
837 nn = index3(n)
838C str utilise com tableau temporaire au lieu de scalaires temporaires
839 str(1,nn)=vec(7,nn)
840 str(2,nn)=vec(8,nn)
841 str(3,nn)=vec(9,nn)
842 ENDDO
843 DO n = 1, nindex3
844 nn = index3(n)
845 vec(7,nn)=vec(1,nn)
846 vec(8,nn)=vec(2,nn)
847 vec(9,nn)=vec(3,nn)
848 vec(1,nn)=-str(1,nn)
849 vec(2,nn)=-str(2,nn)
850 vec(3,nn)=-str(3,nn)
851 ENDDO
852 RETURN
853 END
my_real
#define my_real
Definition cppsort.cpp:32
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dreh
subroutine dreh(b, tr, ii, jj, ken)
Definition matrix.F:161
valpvecdp
subroutine valpvecdp(sig, val, vec, nel)
Definition matrix.F:533
jacobiew
subroutine jacobiew(a, n, ew, ev, nrot)
Definition matrix.F:29
valpvec
subroutine valpvec(sig, val, vec, nel)
Definition matrix.F:215
floatmin
void floatmin(int *a, int *b, float *flm)
Definition precision.c:71
starter
program starter
Definition starter.F:39