CalculiX  2.13
A Free Software Three-Dimensional Structural Finite Element Program
Functions/Subroutines
dlz.f File Reference

Go to the source code of this file.

Functions/Subroutines

subroutine dlzhes (N, A, NA, B, NB, X, NX, WANTX)
 
subroutine dlzit (N, A, NA, B, NB, X, NX, WANTX, ITER, EIGA, EIGB)
 

Function/Subroutine Documentation

◆ dlzhes()

subroutine dlzhes ( integer  N,
complex*16, dimension(na,n)  A,
integer  NA,
complex*16, dimension(nb,n)  B,
integer  NB,
complex*16, dimension(nx,n)  X,
integer  NX,
logical  WANTX 
)
7 C THIS SUBROUTINE REDUCES THE COMPLEX MATRIX A TO UPPER
8 C HESSENBERG FORM AND REDUCES THE COMPLEX MATRIX B TO
9 C TRIANGULAR FORM
10 C INPUT PARAMETERS..
11 C N THE ORDER OF THE A AND B MATRICES
12 C A A COMPLEX MATRIX
13 C NA THE ROW DIMENSION OF THE A MATRIX
14 C B A COMPLEX MATRIX
15 C NB THE ROW DIMENSION OF THE B MATRIX
16 C NX THE ROW DIMENSION OF THE X MATRIX
17 C WANTX A LOGICAL VARIABLE WHICH IS SET TO .TRUE. IF
18 C THE EIGENVECTORS ARE WANTED. OTHERWISE IT SHOULD
19 C BE SET TO .FALSE.
20 C OUTPUT PARAMETERS..
21 C A ON OUTPUT A IS AN UPPER HESSENBERG MATRIX, THE
22 C ORIGINAL MATRIX HAS BEEN DESTROYED
23 C B AN UPPER TRIANGULAR MATRIX, THE ORIGINAL MATRIX
24 C HAS BEEN DESTROYED
25 C X CONTAINS THE TRANSFORMATIONS NEEDED TO COMPUTE
26 C THE EIGENVECTORS OF THE ORIGINAL SYSTEM
27  implicit none
28 !
29  LOGICAL wantx
30 !
31  integer k,nm1,nm2,jm2,jp1,ip1,imj,j,ii,im1,i,nx,nb,na,n
32 !
33  REAL*8 c, d
34 !
35  COMPLEX*16 y, w, z, a(na,n), b(nb,n), x(nx,n)
36 !
37  nm1 = n - 1
38 C REDUCE B TO TRIANGULAR FORM USING ELEMENTARY
39 C TRANSFORMATIONS
40  DO 80 i=1,nm1
41  d = 0.d0
42  ip1 = i + 1
43  DO 10 k=ip1,n
44  y = b(k,i)
45  c = dabs(dreal(y)) + dabs(dimag(y))
46  IF (c.LE.d) GO TO 10
47  d = c
48  ii = k
49  10 CONTINUE
50  IF (d.EQ.0.d0) GO TO 80
51  y = b(i,i)
52  IF (d.LE.dabs(dreal(y))+dabs(dimag(y))) GO TO 40
53 C MUST INTERCHANGE
54  DO 20 j=1,n
55  y = a(i,j)
56  a(i,j) = a(ii,j)
57  a(ii,j) = y
58  20 CONTINUE
59  DO 30 j=i,n
60  y = b(i,j)
61  b(i,j) = b(ii,j)
62  b(ii,j) = y
63  30 CONTINUE
64  40 DO 70 j=ip1,n
65  y = b(j,i)/b(i,i)
66  IF (dreal(y).EQ.0.d0 .AND. dimag(y).EQ.0.d0) GO TO 70
67  DO 50 k=1,n
68  a(j,k) = a(j,k) - y*a(i,k)
69  50 CONTINUE
70  DO 60 k=ip1,n
71  b(j,k) = b(j,k) - y*b(i,k)
72  60 CONTINUE
73  70 CONTINUE
74  b(ip1,i) = (0.d0,0.d0)
75  80 CONTINUE
76 C INITIALIZE X
77  IF (.NOT.wantx) GO TO 110
78  DO 100 i=1,n
79  DO 90 j=1,n
80  x(i,j) = (0.d0,0.d0)
81  90 CONTINUE
82  x(i,i) = (1.d0,0.d0)
83  100 CONTINUE
84 C REDUCE A TO UPPER HESSENBERG FORM
85  110 nm2 = n - 2
86  IF (nm2.LT.1) GO TO 270
87  DO 260 j=1,nm2
88  jm2 = nm1 - j
89  jp1 = j + 1
90  DO 250 ii=1,jm2
91  i = n + 1 - ii
92  im1 = i - 1
93  imj = i - j
94  w = a(i,j)
95  z = a(im1,j)
96  IF (dabs(dreal(w))+dabs(dimag(w)).LE.dabs(dreal(z))
97  * +dabs(dimag(z))) GO TO 140
98 C MUST INTERCHANGE ROWS
99  DO 120 k=j,n
100  y = a(i,k)
101  a(i,k) = a(im1,k)
102  a(im1,k) = y
103  120 CONTINUE
104  DO 130 k=im1,n
105  y = b(i,k)
106  b(i,k) = b(im1,k)
107  b(im1,k) = y
108  130 CONTINUE
109  140 z = a(i,j)
110  IF (dreal(z).EQ.0.d0 .AND. dimag(z).EQ.0.d0) GO TO 170
111  y = z/a(im1,j)
112  DO 150 k=jp1,n
113  a(i,k) = a(i,k) - y*a(im1,k)
114  150 CONTINUE
115  DO 160 k=im1,n
116  b(i,k) = b(i,k) - y*b(im1,k)
117  160 CONTINUE
118 C TRANSFORMATION FROM THE RIGHT
119  170 w = b(i,im1)
120  z = b(i,i)
121  IF (dabs(dreal(w))+dabs(dimag(w)).LE.dabs(dreal(z))
122  * +dabs(dimag(z))) GO TO 210
123 C MUST INTERCHANGE COLUMNS
124  DO 180 k=1,i
125  y = b(k,i)
126  b(k,i) = b(k,im1)
127  b(k,im1) = y
128  180 CONTINUE
129  DO 190 k=1,n
130  y = a(k,i)
131  a(k,i) = a(k,im1)
132  a(k,im1) = y
133  190 CONTINUE
134  IF (.NOT.wantx) GO TO 210
135  DO 200 k=imj,n
136  y = x(k,i)
137  x(k,i) = x(k,im1)
138  x(k,im1) = y
139  200 CONTINUE
140  210 z = b(i,im1)
141  IF (dreal(z).EQ.0.d0 .AND. dimag(z).EQ.0.d0) GO TO 250
142  y = z/b(i,i)
143  DO 220 k=1,im1
144  b(k,im1) = b(k,im1) - y*b(k,i)
145  220 CONTINUE
146  b(i,im1) = (0.d0,0.d0)
147  DO 230 k=1,n
148  a(k,im1) = a(k,im1) - y*a(k,i)
149  230 CONTINUE
150  IF (.NOT.wantx) GO TO 250
151  DO 240 k=imj,n
152  x(k,im1) = x(k,im1) - y*x(k,i)
153  240 CONTINUE
154  250 CONTINUE
155  a(jp1+1,j) = (0.d0,0.d0)
156  260 CONTINUE
157  270 RETURN

◆ dlzit()

subroutine dlzit (   N,
complex*16, dimension(na,n)  A,
  NA,
complex*16, dimension(nb,n)  B,
  NB,
complex*16, dimension(nx,n)  X,
  NX,
logical  WANTX,
integer, dimension(n)  ITER,
complex*16, dimension(n)  EIGA,
complex*16, dimension(n)  EIGB 
)
161 C THIS SUBROUTINE SOLVES THE GENERALIZED EIGENVALUE PROBLEM
162 C A X = LAMBDA B X
163 C WHERE A IS A COMPLEX UPPER HESSENBERG MATRIX OF
164 C ORDER N AND B IS A COMPLEX UPPER TRIANGULAR MATRIX OF ORDER N
165 C INPUT PARAMETERS
166 C N ORDER OF A AND B
167 C A AN N X N UPPER HESSENBERG COMPLEX MATRIX
168 C NA THE ROW DIMENSION OF THE A MATRIX
169 C B AN N X N UPPER TRIANGULAR COMPLEX MATRIX
170 C NB THE ROW DIMENSION OF THE B MATRIX
171 C X CONTAINS TRANSFORMATIONS TO OBTAIN EIGENVECTORS OF
172 C ORIGINAL SYSTEM. IF EIGENVECTORS ARE REQUESTED AND QZHES
173 C IS NOT CALLED, X SHOULD BE SET TO THE IDENTITY MATRIX
174 C NX THE ROW DIMENSION OF THE X MATRIX
175 C WANTX LOGICAL VARIABLE WHICH SHOULD BE SET TO .TRUE.
176 C IF EIGENVECTORS ARE WANTED. OTHERWISE IT
177 C SHOULD BE SET TO .FALSE.
178 C OUTPUT PARAMETERS
179 C X THE ITH COLUMN CONTAINS THE ITH EIGENVECTOR
180 C IF EIGENVECTORS ARE REQUESTED
181 C ITER AN INTEGER ARRAY OF LENGTH N WHOSE ITH ENTRY
182 C CONTAINS THE NUMBER OF ITERATIONS NEEDED TO FIND
183 C THE ITH EIGENVALUE. FOR ANY I IF ITER(I) =-1 THEN
184 C AFTER 30 ITERATIONS THERE HAS NOT BEEN A SUFFICIENT
185 C DECREASE IN THE LAST SUBDIAGONAL ELEMENT OF A
186 C TO CONTINUE ITERATING.
187 C EIGA A COMPLEX ARRAY OF LENGTH N CONTAINING THE DIAGONAL OF A
188 C EIGB A COMPLEX ARRAY OF LENGTH N CONTAINING THE DIAGONAL OF B
189 C THE ITH EIGENVALUE CAN BE FOUND BY DIVIDING EIGA(I) BY
190 C EIGB(I). WATCH OUT FOR EIGB(I) BEING ZERO
191  COMPLEX*16 a(na,n), b(nb,n), eiga(n), eigb(n)
192  COMPLEX*16 s, w, y, z, cdsqrt
193  COMPLEX*16 x(nx,n)
194  INTEGER iter(n)
195  COMPLEX*16 annm1, alfm, betm, d, sl, den, num, anm1m1
196  REAL*8 epsa, epsb, ss, r, anorm, bnorm, ani, bni, c
197  REAL*8 d0, d1, d2, e0, e1
198  LOGICAL wantx
199  nn = n
200 C COMPUTE THE MACHINE PRECISION TIMES THE NORM OF A AND B
201  anorm = 0.d0
202  bnorm = 0.d0
203  DO 30 i=1,n
204  ani = 0.d0
205  IF (i.EQ.1) GO TO 10
206  y = a(i,i-1)
207  ani = ani + dabs(dreal(y)) + dabs(dimag(y))
208  10 bni = 0.d0
209  DO 20 j=i,n
210  ani = ani + dabs(dreal(a(i,j))) + dabs(dimag(a(i,j)))
211  bni = bni + dabs(dreal(b(i,j))) + dabs(dimag(b(i,j)))
212  20 CONTINUE
213  IF (ani.GT.anorm) anorm = ani
214  IF (bni.GT.bnorm) bnorm = bni
215  30 CONTINUE
216  IF (anorm.EQ.0.d0) anorm = 1.d0
217  IF (bnorm.EQ.0.d0) bnorm = 1.d0
218  epsb = bnorm
219  epsa = anorm
220  40 epsa = epsa/2.d0
221  epsb = epsb/2.d0
222  c = anorm + epsa
223  IF (c.GT.anorm) GO TO 40
224  IF (n.LE.1) GO TO 320
225  50 its = 0
226  nm1 = nn - 1
227 C CHECK FOR NEGLIGIBLE SUBDIAGONAL ELEMENTS
228  60 d2 = dabs(dreal(a(nn,nn))) + dabs(dimag(a(nn,nn)))
229  DO 70 lb=2,nn
230  l = nn + 2 - lb
231  ss = d2
232  y = a(l-1,l-1)
233  d2 = dabs(dreal(y)) + dabs(dimag(y))
234  ss = ss + d2
235  y = a(l,l-1)
236  r = ss + dabs(dreal(y)) + dabs(dimag(y))
237  IF (r.EQ.ss) GO TO 80
238  70 CONTINUE
239  l = 1
240  80 IF (l.EQ.nn) GO TO 320
241  IF (its.LT.30) GO TO 90
242  iter(nn) = -1
243  IF (dabs(dreal(a(nn,nm1)))+dabs(dimag(a(nn,nm1))).GT.0.8d0*
244  * dabs(dreal(annm1))+dabs(dimag(annm1))) RETURN
245  90 IF (its.EQ.10 .OR. its.EQ.20) GO TO 110
246 C COMPUTE SHIFT AS EIGENVALUE OF LOWER 2 BY 2
247  annm1 = a(nn,nm1)
248  anm1m1 = a(nm1,nm1)
249  s = a(nn,nn)*b(nm1,nm1) - annm1*b(nm1,nn)
250  w = annm1*b(nn,nn)*(a(nm1,nn)*b(nm1,nm1)-b(nm1,nn)*anm1m1)
251  y = (anm1m1*b(nn,nn)-s)/2.
252  z = cdsqrt(y*y+w)
253  IF (dreal(z).EQ.0.d0 .AND. dimag(z).EQ.0.d0) GO TO 100
254  d0 = dreal(y/z)
255  IF (d0.LT.0.d0) z = -z
256  100 den = (y+z)*b(nm1,nm1)*b(nn,nn)
257  IF (dreal(den).EQ.0.d0 .AND. dimag(den).EQ.0.d0) den =
258  * cmplx(epsa,0.d0)
259  num = (y+z)*s - w
260  GO TO 120
261 C AD-HOC SHIFT
262  110 y = a(nm1,nn-2)
263  num = cmplx(dabs(dreal(annm1))+dabs(dimag(annm1)),dabs(dreal(y))
264  * +dabs(dimag(y)))
265  den = (1.d0,0.d0)
266 C CHECK FOR 2 CONSECUTIVE SMALL SUBDIAGONAL ELEMENTS
267  120 IF (nn.EQ.l+1) GO TO 140
268  d2 = dabs(dreal(a(nm1,nm1))) + dabs(dimag(a(nm1,nm1)))
269  e1 = dabs(dreal(annm1)) + dabs(dimag(annm1))
270  d1 = dabs(dreal(a(nn,nn))) + dabs(dimag(a(nn,nn)))
271  nl = nn - (l+1)
272  DO 130 mb=1,nl
273  m = nn - mb
274  e0 = e1
275  y = a(m,m-1)
276  e1 = dabs(dreal(y)) + dabs(dimag(y))
277  d0 = d1
278  d1 = d2
279  y = a(m-1,m-1)
280  d2 = dabs(dreal(y)) + dabs(dimag(y))
281  y = a(m,m)*den - b(m,m)*num
282  d0 = (d0+d1+d2)*(dabs(dreal(y))+dabs(dimag(y)))
283  e0 = e0*e1*(dabs(dreal(den))+dabs(dimag(den))) + d0
284  IF (e0.EQ.d0) GO TO 150
285  130 CONTINUE
286  140 m = l
287  150 CONTINUE
288  its = its + 1
289  w = a(m,m)*den - b(m,m)*num
290  z = a(m+1,m)*den
291  d1 = dabs(dreal(z)) + dabs(dimag(z))
292  d2 = dabs(dreal(w)) + dabs(dimag(w))
293 C FIND L AND M AND SET A=LAM AND B=LBM
294 C NP1 = N + 1
295  lor1 = l
296  nnorn = nn
297  IF (.NOT.wantx) GO TO 160
298  lor1 = 1
299  nnorn = n
300  160 DO 310 i=m,nm1
301  j = i + 1
302 C FIND ROW TRANSFORMATIONS TO RESTORE A TO
303 C UPPER HESSENBERG FORM. APPLY TRANSFORMATIONS
304 C TO A AND B
305  IF (i.EQ.m) GO TO 170
306  w = a(i,i-1)
307  z = a(j,i-1)
308  d1 = dabs(dreal(z)) + dabs(dimag(z))
309  d2 = dabs(dreal(w)) + dabs(dimag(w))
310  IF (d1.EQ.0.d0) GO TO 60
311  170 IF (d2.GT.d1) GO TO 190
312 C MUST INTERCHANGE ROWS
313  DO 180 k=i,nnorn
314  y = a(i,k)
315  a(i,k) = a(j,k)
316  a(j,k) = y
317  y = b(i,k)
318  b(i,k) = b(j,k)
319  b(j,k) = y
320  180 CONTINUE
321  IF (i.GT.m) a(i,i-1) = a(j,i-1)
322  IF (d2.EQ.0.d0) GO TO 220
323 C THE SCALING OF W AND Z IS DESIGNED TO AVOID A DIVISION BY ZERO
324 C WHEN THE DENOMINATOR IS SMALL
325  y = cmplx(dreal(w)/d1,dimag(w)/d1)/cmplx(dreal(z)/d1,dimag(z)/
326  * d1)
327  GO TO 200
328  190 y = cmplx(dreal(z)/d2,dimag(z)/d2)/cmplx(dreal(w)/d2,dimag(w)/
329  * d2)
330  200 DO 210 k=i,nnorn
331  a(j,k) = a(j,k) - y*a(i,k)
332  b(j,k) = b(j,k) - y*b(i,k)
333  210 CONTINUE
334  220 IF (i.GT.m) a(j,i-1) = (0.d0,0.d0)
335 C PERFORM TRANSFORMATIONS FROM RIGHT TO RESTORE B TO
336 C TRIANGLULAR FORM
337 C APPLY TRANSFORMATIONS TO A
338  z = b(j,i)
339  w = b(j,j)
340  d2 = dabs(dreal(w)) + dabs(dimag(w))
341  d1 = dabs(dreal(z)) + dabs(dimag(z))
342  IF (d1.EQ.0.d0) GO TO 60
343  IF (d2.GT.d1) GO TO 270
344 C MUST INTERCHANGE COLUMNS
345  DO 230 k=lor1,j
346  y = a(k,j)
347  a(k,j) = a(k,i)
348  a(k,i) = y
349  y = b(k,j)
350  b(k,j) = b(k,i)
351  b(k,i) = y
352  230 CONTINUE
353  IF (i.EQ.nm1) GO TO 240
354  y = a(j+1,j)
355  a(j+1,j) = a(j+1,i)
356  a(j+1,i) = y
357  240 IF (.NOT.wantx) GO TO 260
358  DO 250 k=1,n
359  y = x(k,j)
360  x(k,j) = x(k,i)
361  x(k,i) = y
362  250 CONTINUE
363  260 IF (d2.EQ.0.d0) GO TO 310
364  z = cmplx(dreal(w)/d1,dimag(w)/d1)/cmplx(dreal(z)/d1,dimag(z)/
365  * d1)
366  GO TO 280
367  270 z = cmplx(dreal(z)/d2,dimag(z)/d2)/cmplx(dreal(w)/d2,dimag(w)/
368  * d2)
369  280 DO 290 k=lor1,j
370  a(k,i) = a(k,i) - z*a(k,j)
371  b(k,i) = b(k,i) - z*b(k,j)
372  290 CONTINUE
373  b(j,i) = (0.d0,0.d0)
374  IF (i.LT.nm1) a(i+2,i) = a(i+2,i) - z*a(i+2,j)
375  IF (.NOT.wantx) GO TO 310
376  DO 300 k=1,n
377  x(k,i) = x(k,i) - z*x(k,j)
378  300 CONTINUE
379  310 CONTINUE
380  GO TO 60
381  320 CONTINUE
382  eiga(nn) = a(nn,nn)
383  eigb(nn) = b(nn,nn)
384  IF (nn.EQ.1) GO TO 330
385  iter(nn) = its
386  nn = nm1
387  IF (nn.GT.1) GO TO 50
388  iter(1) = 0
389  GO TO 320
390 C FIND EIGENVECTORS USING B FOR INTERMEDIATE STORAGE
391  330 IF (.NOT.wantx) RETURN
392  m = n
393  340 CONTINUE
394  alfm = a(m,m)
395  betm = b(m,m)
396  b(m,m) = (1.d0,0.d0)
397  l = m - 1
398  IF (l.EQ.0) GO TO 370
399  350 CONTINUE
400  l1 = l + 1
401  sl = (0.d0,0.d0)
402  DO 360 j=l1,m
403  sl = sl + (betm*a(l,j)-alfm*b(l,j))*b(j,m)
404  360 CONTINUE
405  y = betm*a(l,l) - alfm*b(l,l)
406  IF (dreal(y).EQ.0.d0 .AND. dimag(y).EQ.0.d0) y =
407  * cmplx((epsa+epsb)/2.d0,0.d0)
408  b(l,m) = -sl/y
409  l = l - 1
410  370 IF (l.GT.0) GO TO 350
411  m = m - 1
412  IF (m.GT.0) GO TO 340
413 C TRANSFORM TO ORIGINAL COORDINATE SYSTEM
414  m = n
415  380 CONTINUE
416  DO 400 i=1,n
417  s = (0.d0,0.d0)
418  DO 390 j=1,m
419  s = s + x(i,j)*b(j,m)
420  390 CONTINUE
421  x(i,m) = s
422  400 CONTINUE
423  m = m - 1
424  IF (m.GT.0) GO TO 380
425 C NORMALIZE SO THAT LARGEST COMPONENT = 1.
426  m = n
427  410 CONTINUE
428  ss = 0.d0
429  DO 420 i=1,n
430  r = dabs(dreal(x(i,m))) + dabs(dimag(x(i,m)))
431  IF (r.LT.ss) GO TO 420
432  ss = r
433  d = x(i,m)
434  420 CONTINUE
435  IF (ss.EQ.0.d0) GO TO 440
436  DO 430 i=1,n
437  x(i,m) = x(i,m)/d
438  430 CONTINUE
439  440 m = m - 1
440  IF (m.GT.0) GO TO 410
441  RETURN
Hosted by OpenAircraft.com, (Michigan UAV, LLC)