V3FIT
blk3d.f
1  SUBROUTINE blk3d(a, bm1, bp1, srces, mblk, nblocks)
2 C-----------------------------------------------
3 C M o d u l e s
4 C-----------------------------------------------
5  USE stel_kinds
6  USE safe_open_mod
7  IMPLICIT NONE
8 C-----------------------------------------------
9 C L o c a l P a r a m e t e r s
10 C-----------------------------------------------
11  INTEGER, PARAMETER :: bytes_per_rprec = 8
12 C-----------------------------------------------
13 C D u m m y A r g u m e n t s
14 C-----------------------------------------------
15  INTEGER, INTENT(in) :: nblocks, mblk
16  REAL(rprec), DIMENSION(mblk,mblk,nblocks), INTENT(in) ::
17  1 a, bp1, bm1
18  REAL(rprec), DIMENSION(mblk,nblocks), INTENT(inout) :: srces
19 C-----------------------------------------------
20 C L o c a l V a r i a b l e s
21 C-----------------------------------------------
22 ! INTEGER :: ibuph, incnow, irecl, incbu, iunit=102, ndisk
23  INTEGER :: k, k1, ier
24  INTEGER, ALLOCATABLE, DIMENSION(:) :: ipiv
25  REAL(rprec), DIMENSION(:,:), ALLOCATABLE :: ainv
26  REAL(rprec), DIMENSION(:,:,:), ALLOCATABLE :: ql
27 C-----------------------------------------------
28 C E x t e r n a l S u b r o u t i n e s
29 C-----------------------------------------------
30  EXTERNAL la_getrf, la_getrs
31 C-----------------------------------------------
32 c modified (June, 2003, ORNL): S. P. Hirshman
33 c-----------------------------------------------------------------------
34 c
35 c this subroutine solves a block-tridiagonal system of equations.
36 c
37 c-----------------------------------------------------------------------
38 c INPUT
39 c mblk : block size
40 c nblocks : number of blocks
41 c a : diagonal block
42 c bm1, bp1 : lower, upper block (see equation below)
43 c srces(mblk,nblocks) : right-hand size source
44 c
45 c LOCAL VARIABLES
46 c iunit : unit number for block-tridiagonal solution disk file.
47 c
48 c solutions are indexed in m-n fourier-space, legendre-space. the tri-diagonal
49 c equation is:
50 c
51 c bm1 * f(l-1) + a * f(l) + bp1 * f(l+1) = source(l)
52 c
53 c GENERAL SOLUTION SCHEME APPLIED TO EACH BLOCK ROW (INDEX L)
54 c
55 c 1. Start from row N and solve for x(N) in terms of x(N-1):
56 c
57 c x(N) = -q(N)*x(N-1) + r(N)
58 c
59 c q(N) = a(N)[-1] * bm1; r(N) = a(N)[-1] * s(N)
60 c
61 c where a(N)[-1] is the inverse of a(N)
62 c
63 c 2. Substitute for lth row to get recursion equation fo q(l) and r(l):
64 c
65 c x(l) = -q(l)*x(l-1) + r(l), in general, where:
66 c
67 c q(l) = (a(l) - bp1(l)*q(l+1))[-1] * bm1(l)
68 c
69 c r(l) = (a(l) - bp1(l)*q(l+1))[-1] * (s(l) - bp1(l)*r(l+1))
70 c
71 c 3. At row l = 1, bm1(1) = 0 and get an equation for x(1) corresponding to q(1) = 0:
72 c
73 c x(1) = r(1)
74 c
75 c 4. Finally, can back-solve for x(N) in terms of x(N-1) from eqn.(1) above
76 c
77 c
78 c NUMERICAL IMPLEMENTATION (USING LAPACK, la_... ROUTINES)
79 c
80 c 1. CALL la_getrf: Perform LU factorization of diagonal block (A)
81 c 2. CALL la_getrs: With multiple (mblk) right-hand sides, to do block inversion
82 c operation, A X = B (store result in X; here B is a matrix)
83 c 3. CALL la_getrs: With single right hand side (source) to solve A x = b (b a vector)
84 c
85 ! ndisk = mblk
86 
87  ALLOCATE (ql(mblk,mblk, nblocks), ainv(mblk,mblk), ipiv(mblk),
88  1 stat=ier)
89  IF (ier .ne. 0) stop 'Allocation error in blk3d'
90 
91 c create disk file for doing direct access i/o.
92 
93 ! incnow = ndisk
94 ! irecl = bytes_per_rprec*incnow
95 ! incbu = 1 + (ndisk - 1)/incnow
96 ! ibuph = 0
97 
98 ! iunit = 10
99 ! CALL safe_open(iunit, ier, 'NULL', 'scratch', 'unformatted',
100 ! 1 irecl, 'DIRECT')
101 ! IF (ier .ne. 0) STOP 'Error opening scratch file in blk3d'
102 
103  ainv = a(:,:,nblocks)
104 
105 c main loop. load and process (backwards) block-rows nblocks to 1.
106 
107  blocks: DO k = nblocks, 1, -1
108 !
109 ! Compute (and save) ql(nblocks) = a(nblocks)[-1] * ql, and source terms A-1 * srces [A-1 == inv(A)]
110 !
111  CALL la_getrf (mblk, mblk, ainv, mblk, ipiv, ier)
112  IF (ier .ne. 0) GOTO 200
113 
114  IF (k .gt. 1) THEN
115  ql(:,:,k) = bm1(:,:,k)
116  CALL la_getrs('n', mblk, mblk, ainv,
117  1 mblk, ipiv, ql(1,1,k), mblk, ier)
118  CALL la_getrs('n', mblk, 1, ainv,
119  1 mblk,ipiv,srces(1,k),mblk,ier)
120 
121 ! CALL wrdisk(iunit, ql, ndisk, incnow, ibuph, incbu, ier)
122 ! IF (ier .ne. 0) GOTO 302
123  ELSE
124  CALL la_getrs('n', mblk, 1, ainv,
125  1 mblk,ipiv,srces(1,k),mblk,ier)
126  END IF
127 
128  k1 = k - 1
129  IF (k1 .eq. 0) EXIT
130 
131 !
132 ! Update effective diagonal "a" matrix and source terms and store ql 2 l-steps back
133 !
134  ainv = a(:,:,k1) - matmul(bp1(:,:,k1), ql(:,:,k))
135  srces(:,k1) = srces(:,k1) - matmul(bp1(:,:,k1),srces(:,k))
136 
137  END DO blocks
138 
139 c backward solution sweep for block-rows k = 2 to nblocks
140 c read blocks ql from disk.
141 
142  DO k = 2, nblocks
143 ! CALL rddisk (iunit, ql, ndisk, incnow, ibuph, ier)
144 ! IF (ier .ne. 0) GOTO 303
145 ! ibuph = ibuph - incbu
146  srces(:,k) = srces(:,k) - matmul(ql(:,:,k),srces(:,k-1))
147  END DO
148 
149  GOTO 400
150 
151 c error returns. ------------------------------------------------------
152 
153  200 CONTINUE
154  WRITE (6, '(a,i4)') ' Error factoring matrix in blk3d: block = '
155  1 , k,' error id = ', ier
156  GOTO 400
157  301 CONTINUE
158 ! WRITE (6, '(a,i8)') ' BLK3D: error in opening file: ',
159 ! 1 'RECL = ', irecl
160  302 CONTINUE
161  WRITE (6, '(a)') ' BLK3D: error in I/O routine WRDISK'
162  303 CONTINUE
163  WRITE (6, '(a)') ' BLK3D: error in I/O routine RDDISK'
164  ier = -2
165  305 CONTINUE
166  WRITE (6, '(2/a,i4,2/)') ' BLK3D: error detected: ier =',
167  1 ier
168  stop
169 
170 c destroy disk file and return. ---------------------------------------
171 
172  400 CONTINUE
173 
174 ! CLOSE (iunit)
175 
176  DEALLOCATE (ainv, ql, ipiv)
177 
178  END SUBROUTINE blk3d
Generated on Thu Mar 5 2020 15:49:23 for V3FIT by   doxygen 1.8.17