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 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,
240 .
241 .
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 MAIN STRESS DEVIATOR
249C . . . . . . . . . . . . . . . . . . .
250 mdemi = -half
251 ttpi = acos(mdemi)
252 ftpi = two*ttpi
253C minimum precision depending on REAL or DOUBLE type
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 Compression accuracy strengthening
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 Rewriting to bypass problem on itanium2 comp 9. + latency = 16
503 DO n = 1, nindex3
504 nn = index3(n)
505C Str uses temporary table instead of temporary scalars
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 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,
558 .
559 .
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 MAIN STRESS DEVIATOR
568C . . . . . . . . . . . . . . . . . . .
569 mdemi = -half
570 ttpi = acos(mdemi)
571 ftpi = two*ttpi
572C minimum precision depending on REAL or DOUBLE type
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 precision reinforcement in simple compression ----
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 Rewriting to bypass problem on itanium2 comp 9. + latency = 16
836 DO n = 1, nindex3
837 nn = index3(n)
838C Str uses temporary table instead of temporary scalars
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
sigeps38
subroutine sigeps38(nel, nuparam, nuvar, nfunc, ifunc, npf, tf, time, timestep, uparam, rho0, rho, volume, eint, epspxx, epspyy, epspzz, epspxy, epspyz, epspzx, depsxx, depsyy, depszz, depsxy, depsyz, depszx, epsxx, epsyy, epszz, epsxy, epsyz, epszx, sigoxx, sigoyy, sigozz, sigoxy, sigoyz, sigozx, signxx, signyy, signzz, signxy, signyz, signzx, sigvxx, sigvyy, sigvzz, sigvxy, sigvyz, sigvzx, soundsp, viscmax, uvar, off, uvarbid, ismstr, mfxx, mfxy, mfxz, mfyx, mfyy, mfyz, mfzx, mfzy, mfzz)
Definition sigeps38.F:48
starter
program starter
Definition starter.F:39