10 real(kind=kreal),
parameter,
private :: pi=3.14159265358979d0
16 integer :: var,dimn,syze,bite
18 var = var + dimn*syze*bite
24 integer,
intent(in) :: n
25 integer,
intent(in),
optional :: n1
26 character(len=*),
intent(inout) :: fname
28 character(len=128) :: tmpname, tmp
30 npos = scan( fname,
'.')
31 nlen = len_trim( fname )
32 if( nlen>128 ) stop
"String too long(>128) in append_int2name"
33 if( n>100000 ) stop
"Integer too big>100000 in append_int2name"
36 write( fname,
'(a,i6)') fname(1:nlen),n
38 write( tmp,
'(i6,a)') n,tmpname(npos:nlen)
39 fname = tmpname(1:npos-1) // adjustl(tmp)
41 if(
present(n1).and.n1/=0)
then
43 fname = fname(1:len_trim(fname))//
'.'//adjustl(tmp)
49 integer,
intent(in) :: iin
50 integer,
pointer :: carray(:)
53 integer,
pointer :: dumarray(:) => null()
54 if( .not.
associated(carray) )
then
58 oldsize =
size( carray )
59 allocate( dumarray(oldsize) )
61 dumarray(i) = carray(i)
64 allocate( carray(oldsize+1) )
66 carray(i) = dumarray(i)
68 carray(oldsize+1) = iin
70 if(
associated(dumarray) )
deallocate( dumarray )
75 real(kind=kreal),
intent(in) :: tensor(6)
76 real(kind=kreal),
intent(out) :: eigval(3)
78 real(kind=kreal) :: i1,i2,i3,r,sita,q, x(3,3), xx(3,3), ii(3,3)
81 ii(1,1)=1.d0; ii(2,2)=1.d0; ii(3,3)=1.d0
82 x(1,1)=tensor(1); x(2,2)=tensor(2); x(3,3)=tensor(3)
83 x(1,2)=tensor(4); x(2,1)=x(1,2)
84 x(2,3)=tensor(5); x(3,2)=x(2,3)
85 x(3,1)=tensor(6); x(1,3)=x(3,1)
88 i1= x(1,1)+x(2,2)+x(3,3)
89 i2= 0.5d0*( i1*i1 - (xx(1,1)+xx(2,2)+xx(3,3)) )
90 i3= x(1,1)*x(2,2)*x(3,3)+x(2,1)*x(3,2)*x(1,3)+x(3,1)*x(1,2)*x(2,3) &
91 -x(3,1)*x(2,2)*x(1,3)-x(2,1)*x(1,2)*x(3,3)-x(1,1)*x(3,2)*x(2,3)
93 r=(-2.d0*i1*i1*i1+9.d0*i1*i2-27.d0*i3)/54.d0
94 q=(i1*i1-3.d0*i2)/9.d0
95 sita = acos(r/dsqrt(q*q*q))
97 eigval(1) = -2.d0*q*cos(sita/3.d0)+i1/3.d0
98 eigval(2) = -2.d0*q*cos((sita+2.d0*pi)/3.d0)+i1/3.d0
99 eigval(3) = -2.d0*q*cos((sita-2.d0*pi)/3.d0)+i1/3.d0
105 subroutine eigen3 (tensor, eigval, princ)
106 real(kind=kreal),
intent(in) :: tensor(6)
107 real(kind=kreal),
intent(out) :: eigval(3)
108 real(kind=kreal),
intent(out) :: princ(3, 3)
110 integer,
parameter :: msweep = 50
111 integer :: i,j, is, ip, iq, ir
112 real(kind=kreal) :: fsum, od, theta, t, c, s, tau, g, h, hd, btens(3,3), factor
114 btens(1,1)=tensor(1); btens(2,2)=tensor(2); btens(3,3)=tensor(3)
115 btens(1,2)=tensor(4); btens(2,1)=btens(1,2)
116 btens(2,3)=tensor(5); btens(3,2)=btens(2,3)
117 btens(3,1)=tensor(6); btens(1,3)=btens(3,1)
127 eigval(i) = btens(i, i)
128 factor = factor + dabs(btens(i,i))
131 if( factor == 0.d0 .or. factor /= factor )
then
134 eigval(1:3) = eigval(1:3)/factor
135 btens(1:3,1:3) = btens(1:3,1:3)/factor
145 fsum = fsum + abs( btens(ip, iq) )
151 if ( fsum < 1.d-10 )
then
152 eigval(1:3) = eigval(1:3)*factor
160 od = 100.d0 * abs(btens(ip, iq) )
161 if ( (od+abs(eigval(ip) ) /= abs(eigval(ip) )) &
162 .and. (od+abs(eigval(iq) ) /= abs(eigval(iq) )))
then
163 hd = eigval(iq) - eigval(ip)
167 if ( abs(hd) + od == abs(hd) )
then
168 t = btens(ip, iq) / hd
170 theta = 0.5d0 * hd / btens(ip, iq)
171 t = 1.d0 / (abs(theta) + sqrt(1.d0 + theta**2) )
172 if ( theta < 0.d0 ) t = - t
177 c = 1.d0 / sqrt(1.d0 + t**2)
180 h = t * btens(ip, iq)
181 eigval(ip) = eigval(ip) - h
182 eigval(iq) = eigval(iq) + h
187 g = btens(min(ir, ip), max(ir, ip) )
188 h = btens(min(ir, iq), max(ir, iq) )
189 btens(min(ir, ip), max(ir, ip) ) = g &
191 btens(min(ir, iq), max(ir, iq) ) = h &
199 princ(ir, ip) = g - s * (h + g * tau)
200 princ(ir, iq) = h + s * (g - h * tau)
210 stop
' Jacobi iteration unable to converge'
215 real(kind=kreal) :: mat(6)
216 real(kind=kreal) :: xj(3,3)
218 xj(1,1)=mat(1); xj(2,2)=mat(2); xj(3,3)=mat(3)
219 xj(1,2)=mat(4); xj(2,1)=xj(1,2)
220 xj(2,3)=mat(5); xj(3,2)=xj(2,3)
221 xj(3,1)=mat(6); xj(1,3)=xj(3,1)
224 +xj(2,1)*xj(3,2)*xj(1,3) &
225 +xj(3,1)*xj(1,2)*xj(2,3) &
226 -xj(3,1)*xj(2,2)*xj(1,3) &
227 -xj(2,1)*xj(1,2)*xj(3,3) &
228 -xj(1,1)*xj(3,2)*xj(2,3)
233 real(kind=kreal) :: xj(3,3)
236 +xj(2,1)*xj(3,2)*xj(1,3) &
237 +xj(3,1)*xj(1,2)*xj(2,3) &
238 -xj(3,1)*xj(2,2)*xj(1,3) &
239 -xj(2,1)*xj(1,2)*xj(3,3) &
240 -xj(1,1)*xj(3,2)*xj(2,3)
245 character(*) :: sub_name
246 integer(kind=kint) :: imsg
247 integer(kind=kint) :: ierr
250 write(imsg,*)
'Memory overflow at ', sub_name
251 write(*,*)
'Memory overflow at ', sub_name
252 call hecmw_abort( hecmw_comm_get_comm( ) )
258 integer,
intent(in) :: NN
259 real(kind=kreal),
intent(inout) :: a(nn,nn)
261 integer :: I, J,K,IW,LR,IP(NN)
262 real(kind=kreal) :: w,wmax,pivot,api,eps,det
281 write(*,
'(''PIVOT ERROR AT'',I5)') k
304 if (j.NE.k) a(i,j)=a(i,j)-w*a(k,j)
330 real(kind=kreal),
intent(in) :: v1(3),v2(3)
331 real(kind=kreal),
intent(out) :: vn(3)
333 vn(1) = v1(2)*v2(3) - v1(3)*v2(2)
334 vn(2) = v1(3)*v2(1) - v1(1)*v2(3)
335 vn(3) = v1(1)*v2(2) - v1(2)*v2(1)
339 real(kind=kreal),
intent(in) :: jacob(3,3)
340 real(kind=kreal),
intent(out) :: tm(6,6)
346 tm(i,j)= jacob(i,j)*jacob(i,j)
348 tm(i,4) = jacob(i,1)*jacob(i,2)
349 tm(i,5) = jacob(i,2)*jacob(i,3)
350 tm(i,6) = jacob(i,3)*jacob(i,1)
352 tm(4,1) = 2.d0*jacob(1,1)*jacob(2,1)
353 tm(5,1) = 2.d0*jacob(2,1)*jacob(3,1)
354 tm(6,1) = 2.d0*jacob(3,1)*jacob(1,1)
355 tm(4,2) = 2.d0*jacob(1,2)*jacob(2,2)
356 tm(5,2) = 2.d0*jacob(2,2)*jacob(3,2)
357 tm(6,2) = 2.d0*jacob(3,2)*jacob(1,2)
358 tm(4,3) = 2.d0*jacob(1,3)*jacob(2,3)
359 tm(5,3) = 2.d0*jacob(2,3)*jacob(3,3)
360 tm(6,3) = 2.d0*jacob(3,3)*jacob(1,3)
361 tm(4,4) = jacob(1,1)*jacob(2,2) + jacob(1,2)*jacob(2,1)
362 tm(5,4) = jacob(2,1)*jacob(3,2) + jacob(2,2)*jacob(3,1)
363 tm(6,4) = jacob(3,1)*jacob(1,2) + jacob(3,2)*jacob(1,1)
364 tm(4,5) = jacob(1,2)*jacob(2,3) + jacob(1,3)*jacob(2,2)
365 tm(5,5) = jacob(2,2)*jacob(3,3) + jacob(2,3)*jacob(3,2)
366 tm(6,5) = jacob(3,2)*jacob(1,3) + jacob(3,3)*jacob(1,2)
367 tm(4,6) = jacob(1,3)*jacob(2,1) + jacob(1,1)*jacob(2,3)
368 tm(5,6) = jacob(2,3)*jacob(3,1) + jacob(2,1)*jacob(3,3)
369 tm(6,6) = jacob(3,3)*jacob(1,1) + jacob(3,1)*jacob(1,3)
377 real(kind=kreal) :: tensor(1:6)
378 real(kind=kreal) :: eigval(3)
379 real(kind=kreal) :: princmatrix(3,3)
380 real(kind=kreal) :: princnormal(3,3)
381 real(kind=kreal) :: tempv(3)
382 real(kind=kreal) :: temps
384 call eigen3(tensor,eigval,princnormal)
386 if (eigval(1)<eigval(2))
then
390 tempv(:)=princnormal(:,1)
391 princnormal(:,1)=princnormal(:,2)
392 princnormal(:,2)=tempv(:)
394 if (eigval(1)<eigval(3))
then
398 tempv(:)=princnormal(:,1)
399 princnormal(:,1)=princnormal(:,3)
400 princnormal(:,3)=tempv(:)
402 if (eigval(2)<eigval(3))
then
406 tempv(:)=princnormal(:,2)
407 princnormal(:,2)=princnormal(:,3)
408 princnormal(:,3)=tempv(:)
413 princmatrix(i,j) = princnormal(i,j) * eigval(j)
422 real(kind=kreal) :: tensor(6)
423 real(kind=kreal) :: eigval(3)
424 real(kind=kreal) :: princ(3,3)
426 real(kind=kreal) :: s11, s22, s33, s12, s23, s13, j1, j2, j3
427 real(kind=kreal) :: ml,nl
428 complex(kind=kreal):: x1,x2,x3
429 real(kind=kreal):: rtemp
440 j2 = -s11*s22 - s22*s33 - s33*s11 + s12**2 + s23**2 + s13**2
441 j3 = s11*s22*s33 + 2*s12*s23*s13 - s11*s23**2 - s22*s13**2 - s33*s12**2
445 call cardano(-j1, -j2, -j3, x1, x2, x3)
449 if (eigval(1)<eigval(2))
then
454 if (eigval(1)<eigval(3))
then
459 if (eigval(2)<eigval(3))
then
466 if (eigval(i)/(eigval(1)+eigval(2)+eigval(3)) < 1.0d-10 )
then
471 ml = ( s23*s13 - s12*(s33-eigval(i)) ) / ( -s23**2 + (s22-eigval(i))*(s33-eigval(i)) )
472 nl = ( s12**2 - (s22-eigval(i))*(s11-eigval(i)) ) / ( s12*s23 - s13*(s22-eigval(i)) )
473 if (abs(ml) >= huge(ml))
then
476 if (abs(nl) >= huge(nl))
then
479 princ(i,1) = eigval(i)/sqrt( 1 + ml**2 + nl**2)
480 princ(i,2) = ml * princ(i,1)
481 princ(i,3) = nl * princ(i,1)
488 real(kind=kreal):: a,b,c
489 real(kind=kreal):: p,q,d
490 complex(kind=kreal):: w
491 complex(kind=kreal):: u,v,y
492 complex(kind=kreal):: x1,x2,x3
493 w = (-1.0d0 + sqrt(dcmplx(-3.0d0)))/2.0d0
494 p = -a**2/9.0d0 + b/3.0d0
495 q = 2.0d0/2.7d1*a**3 - a*b/3.0d0 + c
496 d = q**2 + 4.0d0*p**3
498 u = ((-dcmplx(q) + sqrt(dcmplx(d)))/2.0d0)**(1.0d0/3.0d0)
502 x1 = u + v -dcmplx(a)/3.0d0
503 x2 = u*w + v*w**2 -dcmplx(a)/3.0d0
504 x3 = u*w**2 + v*w -dcmplx(a)/3.0d0
506 y = (-dcmplx(q))**(1.0d0/3.0d0)
507 x1 = y -dcmplx(a)/3.0d0
508 x2 = y*w -dcmplx(a)/3.0d0
509 x3 = y*w**2 -dcmplx(a)/3.0d0
515 real(kind=kreal),
intent(in) :: r(3)
516 real(kind=kreal),
intent(inout) :: v(3)
518 real(kind=kreal) :: rotv(3), rv
519 real(kind=kreal) :: cosx, sinc(2)
520 real(kind=kreal) :: x, x2, x4, x6
521 real(kind=kreal),
parameter :: c0 = 0.5d0
522 real(kind=kreal),
parameter :: c2 = -4.166666666666666d-002
523 real(kind=kreal),
parameter :: c4 = 1.388888888888889d-003
524 real(kind=kreal),
parameter :: c6 = -2.480158730158730d-005
526 x2 = dot_product(r,r)
532 sinc(1) = 1.d0-x2/6.d0+x4/120.d0
533 sinc(2) = c0+c2*x2+c4*x4+c6*x6
536 sinc(2) = (1.d0-cosx)/x2
540 rv = dot_product(r,v)
541 rotv(1:3) = cosx*v(1:3)
542 rotv(1:3) = rotv(1:3)+rv*sinc(2)*r(1:3)
543 rotv(1) = rotv(1) + (-v(2)*r(3)+v(3)*r(2))*sinc(1)
544 rotv(2) = rotv(2) + (-v(3)*r(1)+v(1)*r(3))*sinc(1)
545 rotv(3) = rotv(3) + (-v(1)*r(2)+v(2)*r(1))*sinc(1)
This module provides aux functions.
subroutine eigen3(tensor, eigval, princ)
Compute eigenvalue and eigenvetor for symmetric 3*3 tensor using Jacobi iteration adapted from numeri...
subroutine cross_product(v1, v2, vn)
subroutine eigen3d(tensor, eigval, princ)
real(kind=kreal) function determinant33(XJ)
Compute determinant for 3*3 matrix.
subroutine insert_int2array(iin, carray)
Insert an integer into a integer array.
real(kind=kreal) function determinant(mat)
Compute determinant for symmetric 3*3 matrix.
subroutine append_int2name(n, fname, n1)
Insert an integer at end of a file name.
subroutine tensor_eigen3(tensor, eigval)
Given symmetric 3x3 matrix M, compute the eigenvalues.
subroutine get_principal(tensor, eigval, princmatrix)
subroutine transformation(jacob, tm)
subroutine fstr_chk_alloc(imsg, sub_name, ierr)
subroutine cardano(a, b, c, x1, x2, x3)
subroutine memget(var, dimn, syze)
Record used memeory.
subroutine rotate_3dvector_by_rodrigues_formula(r, v)
subroutine calinverse(NN, A)
calculate inverse of matrix a