1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
|
.help svdsol Jun99 "Slalib Package"
.nf
SUBROUTINE slSVDS (M, N, MP, NP, B, U, W, V, WORK, X)
- - - - - - -
S V D S
- - - - - - -
From a given vector and the SVD of a matrix (as obtained from
the SVD routine), obtain the solution vector (double precision)
This routine solves the equation:
A . x = b
where:
A is a given M (rows) x N (columns) matrix, where M.GE.N
x is the N-vector we wish to find
b is a given M-vector
by means of the Singular Value Decomposition method (SVD). In
this method, the matrix A is first factorised (for example by
the routine slSVD) into the following components:
A = U x W x VT
where:
A is the M (rows) x N (columns) matrix
U is an M x N column-orthogonal matrix
W is an N x N diagonal matrix with W(I,I).GE.0
VT is the transpose of an NxN orthogonal matrix
Note that M and N, above, are the LOGICAL dimensions of the
matrices and vectors concerned, which can be located in
arrays of larger PHYSICAL dimensions MP and NP.
The solution is found from the expression:
x = V . [diag(1/Wj)] . (transpose(U) . b)
Notes:
1) If matrix A is square, and if the diagonal matrix W is not
adjusted, the method is equivalent to conventional solution
of simultaneous equations.
2) If M>N, the result is a least-squares fit.
3) If the solution is poorly determined, this shows up in the
SVD factorisation as very small or zero Wj values. Where
a Wj value is small but non-zero it can be set to zero to
avoid ill effects. The present routine detects such zero
Wj values and produces a sensible solution, with highly
correlated terms kept under control rather than being allowed
to elope to infinity, and with meaningful values for the
other terms.
Given:
M,N i numbers of rows and columns in matrix A
MP,NP i physical dimensions of array containing matrix A
B d(M) known vector b
U d(MP,NP) array containing MxN matrix U
W d(N) NxN diagonal matrix W (diagonal elements only)
V d(NP,NP) array containing NxN orthogonal matrix V
Returned:
WORK d(N) workspace
X d(N) unknown vector x
Reference:
Numerical Recipes, section 2.9.
P.T.Wallace Starlink 29 October 1993
Copyright (C) 1995 Rutherford Appleton Laboratory
Copyright (C) 1995 Association of Universities for Research in Astronomy Inc.
.fi
.endhelp
|
