1 : C PDIRE3 SOURCE CHAT 05/01/13 02:11:26 5004 2 : c--------------------------------------------------------------------- 3 : c 4 : SUBROUTINE PDIRE3 (S, XL, V) 5 : c 6 : c===================================================================== 7 : c = 8 : c This routine computes the eigenvectors. = 9 : c = 10 : c Input : xl (3) eigenvalues = 11 : c s (6) original matrix = 12 : c Output: v (3,3) eigenvectors (normalised) = 13 : c = 14 : c Note: s = (Sxx, Syy, Szz, Sxy, Sxz, Syz) = 15 : c xl = (S11, S22, S33) = 16 : c = 17 : c===================================================================== 18 : IMPLICIT INTEGER(I-N) 19 : real*8 s (6), xl (3), v (3,3) 20 : c 21 : real*8 toler, comp, compi, comp2, a1, a2, a3, dot, xn 22 : c 23 : parameter (toler = 1.0 d-07) 24 : c 25 : comp2 = xl (1)**2 + xl (2)**2 + xl (3)**2 26 : compi = 1.0 d0 / sqrt(comp2) 27 : c 28 : a1 = abs ((xl(2)-xl(3))*compi) 29 : a2 = abs ((xl(1)-xl(3))*compi) 30 : a3 = abs ((xl(1)-xl(2))*compi) 31 : c=====---------------------------------------------------------------- 32 : c case-a: three equal eigenvalues | 33 : c=====---------------------------------------------------------------- 34 : if (a1 .lt. toler .and. a2 .lt. toler) then 35 : v (1,1) = 1.0 d0 36 : v (2,2) = 1.0 d0 37 : v (3,3) = 1.0 d0 38 : v (1,2) = 0.0 d0 39 : v (1,3) = 0.0 d0 40 : v (2,1) = 0.0 d0 41 : v (2,3) = 0.0 d0 42 : v (3,1) = 0.0 d0 43 : v (3,2) = 0.0 d0 44 : return 45 : end if 46 : c=====---------------------------------------------------------------- 47 : c case-b: two equal eigenvalues | 48 : c=====---------------------------------------------------------------- 49 : if (a1 .lt. toler) then 49 : call pvecto (s, xl(1), v(1,1), comp2) 49 : call vecbas (v(1,1), v(1,2), v(1,3)) 50 : return 51 : end if 52 : if (a2 .lt. toler) then 52 : call pvecto (s, xl(2), v(1,2), comp2) 52 : call vecbas (v(1,2), v(1,3), v(1,1)) 53 : return 54 : end if 55 : if (a3 .lt. toler) then 55 : call pvecto (s, xl (3), v(1,3), comp2) 55 : call vecbas (v(1,3), v(1,1), v(1,2)) 56 : return 57 : end if 58 : c=====---------------------------------------------------------------- 59 : c case-c: three different eigenvalues | 60 : c=====---------------------------------------------------------------- 60 : call pvecto (s, xl(1), v(1,1), comp2) 60 : call pvecto (s, xl(2), v(1,2), comp2) 61 : c 62 : dot = v (1,1) * v (1,2) + v (2,1) * v (2,2) + v (3,1) * v (3,2) 63 : c 64 : v (1,2) = v (1,2) - dot * v (1,1) 65 : v (2,2) = v (2,2) - dot * v (2,1) 66 : v (3,2) = v (3,2) - dot * v (3,1) 67 : c 68 : xn = 1.0 d0 / SQRT (v(1,2)**2+v(2,2)**2+v(3,2)**2) 69 : c 70 : v (1,2) = v (1,2) * xn 71 : v (2,2) = v (2,2) * xn 72 : v (3,2) = v (3,2) * xn 73 : v (1,3) = v (2,1) * v (3,2) - v (3,1) * v (2,2) 74 : v (2,3) = v (3,1) * v (1,2) - v (1,1) * v (3,2) 75 : v (3,3) = v (1,1) * v (2,2) - v (2,1) * v (1,2) 76 : c 77 : return 78 : end 79 : 80 :
© Cast3M 2003 - All rights reserved.
Disclaimer