V3FIT
svdinv.f
1  SUBROUTINE svdinv(m, a, wcut, file)
2  USE stel_kinds
3  IMPLICIT NONE
4 C-----------------------------------------------
5 C D u m m y A r g u m e n t s
6 C-----------------------------------------------
7  INTEGER :: m
8  REAL(rprec) :: wcut
9  CHARACTER(LEN=*) :: file
10  REAL(rprec), DIMENSION(m,m) :: a
11 C-----------------------------------------------
12 C L o c a l P a r a m e t e r s
13 C-----------------------------------------------
14  REAL(rprec) :: zero = 0
15 C-----------------------------------------------
16 C L o c a l V a r i a b l e s
17 C-----------------------------------------------
18  INTEGER :: nonzero, nlast, istat, i, j, iunit
19  REAL(rprec), ALLOCATABLE :: u(:,:), w(:), v(:,:)
20  REAL(rprec) :: wmax, wcuta
21 C-----------------------------------------------
22 
23 c Inverts Square Matrix A using SVD method
24 c Allows supression of weights, storage and recall of answers
25 c Uses svd routines from numerical recipes
26 c Prashant Valanju (Sept 1998) pvalanju@mail.utexas.edu
27 
28 c A(M,M) - Matrix on input, Inverse on output
29 c M - Size of A
30 c Wcut - Cutoff parameter
31 c If Wcut = 0, writes answers into "file", all weights are kept
32 c If Wcut < 0, reads previous answers from "file",
33 c and uses -Wcut for weights to keep
34 c If 0 < Wcut < 1 : Cuts off weights with w(i)/wmax < wcut
35 c Use this to cut off small weights, does calculation
36 c If Wcut > 1, keeps INT(Wcut) weights, does calculation
37 c Use this to cut off fixed number of weights
38 c
39 c Note: M is same as N, Mp, and Np in this routine
40 
41 c ALLOCATE local arrays
42  ALLOCATE (u(m,m), w(m), v(m,m), stat=istat)
43  IF (istat .eq. 0) THEN !have enough memory
44 
45  iunit = 41
46  wcuta = abs(wcut)
47 c...............................
48 c If Wcut < 0 use precalculated weights to get answer fast
49 c use wcuta = ABS(wcut) for everything
50 c If any error happens while reading "file", do full svd calculation
51  IF (wcut .lt. zero) THEN !Read file from previous calculation
52  CALL safe_open(iunit, istat, file, 'old', 'formatted')
53  IF (istat .ne. 0) GOTO 98
54  READ (iunit, *, err=98)
55  READ (iunit, *, err=98)
56  READ (iunit, *, err=98) m, nonzero
57 c !Read exactly the same way they were written
58  DO i = 1, nonzero
59  READ (iunit, *, err=98) w(i) !read i-th weight
60  DO j = 1, m
61  READ (iunit, *, err=98) v(j,i), u(j,i)
62  END DO
63  END DO
64  CLOSE(unit=iunit)
65 
66 c Decide how many weights to keep
67  nlast = wcuta !remember Wcut is < 0
68  nlast = min(nonzero,nlast) !Do not use nonzero weights
69 
70  GOTO 99 !bypass svdcmp and GOTO matrix inversion
71  ENDIF !End of precalculated U,w,V branch Wcut < 0
72  98 CONTINUE
73 c.......................................
74 c If Wcut .ge. 0, do full svd calculation and save V,w,U
75 c Initialize all to zero to wipe out effects of old call
76  DO j = 1, m
77  w(j) = zero !Zero all weights
78  u(:m,j) = a(:m,j) !Because U will be changed by svdcmp
79  v(:m,j) = zero
80  END DO
81 
82  CALL svdcmp (u, m, m, m, m, w, v) !Do SVD decomposition
83 
84 c Sort weights and matrices with DECREASING weights
85 c Permute weight w(i) AND column vectors U(*,i), V(*,i) at the same time
86  CALL sortsvd (m, m, m, m, w, u, v)
87 
88 c Find the number of nonzero weights (already decreasing ordered)
89  DO nonzero = m, 1, -1
90 c !Found first nonzero weight, get out
91  IF (w(nonzero) .ne. 0) EXIT
92  END DO
93  IF (nonzero .le. 0) GOTO 999
94 
95 c Write all weights and U, V into SVD file (if Wcut=0)
96  IF (wcut .eq. zero) THEN !no nonzero weights, you must be kidding
97  CALL safe_open(iunit, istat, file, 'unknown', 'formatted')
98  WRITE (iunit, *) 'Max w = ', w(1), ', Min w = ', w(nonzero)
99  WRITE (iunit, *) 'Ratio = ', w(1)/w(nonzero)
100  WRITE (iunit, *) m, nonzero
101 c !WRITE exactly the same way they were written
102  DO i = 1, nonzero
103  WRITE (iunit, *) w(i) !write i-th weight
104  DO j = 1, m
105  WRITE (iunit, *) v(j,i), u(j,i)
106  END DO
107  END DO
108  CLOSE(unit=iunit)
109  END IF
110 c.......................................
111 
112 c Following gets done in both cases (using saved or recalculating)
113 c Decide how many weights to keep
114  99 CONTINUE
115  IF (wcuta .gt. 1) THEN
116  nlast = wcuta
117  ELSE !cutoff small weights
118  wmax = w(1) !weights are already ordered
119  nlast = 0
120  DO i = 1, nonzero
121  IF (w(i) .gt. wmax*wcuta) THEN
122  nlast = nlast + 1 !accept this weight
123  ELSE
124  GOTO 96
125  END IF
126  END DO
127  END IF
128 
129 c Find Inverse of A = V [diag(1/wi)] Utranspose, return in A
130 c First do [diag(1/wi)] Utranspose, store answer in U again
131  96 CONTINUE
132  IF (nlast.le.0 .or. nlast.gt.nonzero) nlast = nonzero
133  DO i = 1, nlast
134 c !divide ith row of Utr by w(i)
135  u(:m,i) = u(:m,i)/w(i)
136  END DO
137 c Zero the infinite 1/weights
138  DO i = nlast + 1, m
139  u(:m,i) = zero !multiply ith row of Utr by zero
140  END DO
141 c Next multiply by V matrix to get A inverse, put it in A
142  DO i = 1, m !multiply ith row vector of V by
143  DO j = 1, m !jth column vector of wU
144  !vector multiply V and wU
145  a(i,j) = sum(v(i,:m)*u(j,:m)) !Put inverse back in A
146  END DO
147  END DO
148 c................................................
149  END IF
150  999 CONTINUE
151 
152  DEALLOCATE (u, w, v, stat=istat)
153 
154  END SUBROUTINE svdinv