mirror of
https://github.com/eddyem/BTA_lib.git
synced 2025-12-06 02:35:20 +03:00
127 lines
3.8 KiB
Fortran
127 lines
3.8 KiB
Fortran
SUBROUTINE sla_SVDSOL (M, N, MP, NP, B, U, W, V, WORK, X)
|
|
*+
|
|
* - - - - - - -
|
|
* S V D S O L
|
|
* - - - - - - -
|
|
*
|
|
* 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 sla_SVD) 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
|
|
*
|
|
* License:
|
|
* This program is free software; you can redistribute it and/or modify
|
|
* it under the terms of the GNU General Public License as published by
|
|
* the Free Software Foundation; either version 2 of the License, or
|
|
* (at your option) any later version.
|
|
*
|
|
* This program is distributed in the hope that it will be useful,
|
|
* but WITHOUT ANY WARRANTY; without even the implied warranty of
|
|
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
|
* GNU General Public License for more details.
|
|
*
|
|
* You should have received a copy of the GNU General Public License
|
|
* along with this program (see SLA_CONDITIONS); if not, write to the
|
|
* Free Software Foundation, Inc., 59 Temple Place, Suite 330,
|
|
* Boston, MA 02111-1307 USA
|
|
*
|
|
*-
|
|
|
|
IMPLICIT NONE
|
|
|
|
INTEGER M,N,MP,NP
|
|
DOUBLE PRECISION B(M),U(MP,NP),W(N),V(NP,NP),WORK(N),X(N)
|
|
|
|
INTEGER J,I,JJ
|
|
DOUBLE PRECISION S
|
|
|
|
|
|
|
|
* Calculate [diag(1/Wj)] . transpose(U) . b (or zero for zero Wj)
|
|
DO J=1,N
|
|
S=0D0
|
|
IF (W(J).NE.0D0) THEN
|
|
DO I=1,M
|
|
S=S+U(I,J)*B(I)
|
|
END DO
|
|
S=S/W(J)
|
|
END IF
|
|
WORK(J)=S
|
|
END DO
|
|
|
|
* Multiply by matrix V to get result
|
|
DO J=1,N
|
|
S=0D0
|
|
DO JJ=1,N
|
|
S=S+V(J,JJ)*WORK(JJ)
|
|
END DO
|
|
X(J)=S
|
|
END DO
|
|
|
|
END
|