V3FIT
svd_solve.f
1  SUBROUTINE svd_solve(m, n, mp, np, a, b, v, w, nw, small)
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, n, mp, np, nw
8  REAL(rprec), DIMENSION(mp,np) :: a
9  REAL(rprec), DIMENSION(n,n) :: v
10  REAL(rprec), DIMENSION(n) :: w
11  REAL(rprec), DIMENSION(m) :: b
12  REAL(rprec) :: small
13 C-----------------------------------------------
14 C L o c a l P a r a m e t e r s
15 C-----------------------------------------------
16 ! REAL(rprec) :: small = 1.0e-10_dp
17 C-----------------------------------------------
18 C L o c a l V a r i a b l e s
19 C-----------------------------------------------
20  INTEGER :: istat, i, j
21  REAL(rprec), ALLOCATABLE :: u(:,:)
22 C-----------------------------------------------
23 
24 c Solves Matrix equation A * V = B for V using SVD method
25 c Uses svd routines from numerical recipes
26 c Prashant Valanju (Dec 1998) pvalanju@mail.utexas.edu
27 
28 c Almost same as SvdSolveB except this one finds
29 c nw solutions V(:,i) by keeping i = 1 to nw vectors, where
30 c nw is the number of vectors with nonzero weights > small
31 
32 c It also uses V internally for SVD decomposition to save space, so
33 
34 c Inputs:
35 c A(M,N) - Matrix with physical dimensions (Mp,Np)
36 c B(M) - R.H.S. of A X = B
37 c M - Number of rows of A to use
38 c N - Number of columns of A to use
39 c Mp,Np - Storage dimensions of A(Mp,Np)
40 c
41 c Output:
42 c nw - number of vectors with normalized weights above small
43 c Eventually, this will be the number of optimum weights
44 c after we decide on a criterion for optimization
45 c V(N,nw) - nw Solutions of A V = B, each a vector of length N
46 c V(:N,iw) = solution with top iw weights kept
47 c Physical dimensions of V are (N,N)
48 c w(n) - Weights in decreasing order
49 c
50 
51 c Allocate local arrays.
52 c Note u(m,n) is enough because needed part of a(mp,np) is copied to
53 c u(m,n), and a(mp,np) is never used directly. This saves space.
54 c It is essential to USE a local u since svdcmp changes it
55 
56  ALLOCATE (u(m,n), stat=istat)
57  IF(istat.ne.0) stop 'Stop: No memory in svd_nesc'
58 
59 c.......................................
60 c Initialize ALL to zero to wipe out effects of old CALL
61  w(:n) = 0 !Zero ALL weights
62  DO j = 1, n
63  u(:m,j) = a(:m,j) !Because U will be changed by svdcmp
64  v(:n,j) = 0
65  END DO
66 
67 c Do the SVD decomposition of a, i.e, of u into u, v and w
68  CALL svdcmp (u, m, n, m, n, w, v)
69 
70 c Sort weights and matrices with DECREASING weights so w(1) is biggest
71 c Permute weight w(i) AND column vectors U(*,i), V(*,i) at the same time
72  CALL sortsvd (m, n, m, n, w, u, v)
73 
74 c Find nw = number of large weights (dcreasing ordered by sortsvd)
75  DO nw = n, 1, -1 !Find first large weight and get out
76  IF ( abs(w(nw)/w(1)) .gt. small) EXIT
77  END DO
78 
79 c Find nw solutions by successively adding up eigenvectors V
80 c with correct weight given by (U(i) dot b)/w(i) in eq 14.3.17 (NR).
81 c The coeff of ith vector V(i) is the dot product of the i-th column
82 c of U with the rhs vector B, divided by w(i)
83 c This does the svdbksb directly, faster than the NR svdbksub routine
84 c and uses less memory due to the dual role of V
85 c Note: any optimization scheme to find the 'best' nw will
86 
87 c First set the 1st vector with largest weight w(1)
88  v(:n,1) = sum(u(:m,1)*b(:m)) *v(:n,1) /w(1)
89 c Next add the vectors with successive weights (in decreasing order)
90  DO i = 2, nw
91  j = i - 1
92  v(:n,i) = v(:n,j) + sum(u(:m,i)*b(:m)) *v(:n,i) / w(i)
93  END DO
94 c................................................
95 
96  DEALLOCATE (u, stat=istat)
97 
98  END SUBROUTINE svd_solve