Upload New File

; $Id: svd_inv_matrix.pro#0.1 $

;

;+

; NAME:

; SVD_INV_MATRIX

;

; PURPOSE:

; Invert a matrix by SVD.

;

; CATEGORY:

; Mathematics.

;

; CALLING SEQUENCE:

; Result = SVD_INV_MATRIX(A, Sv, Tr)

;

; INPUTS:

; A: Input matrix (floating-point, single or double precision).

;

; KEYWORD PARAMETERS:

; TOL: Relative tolerance on singular values of matrix A. 

;    Singular values lower than 

;    tol * MAX(ABS(sv)) are set to 0 in the matrix inversion.

;    Default TOL = 1e-6.

;    

; DOUBLE: set this binary keyword to perform inversion of matrix A 

;    in double-precision. This is the default if A is double-precision.

;

; OUTPUTS:

; Returns inverse of A.

;

; OPTIONAL OUTPUTS:

; Sv: vector of singular values of (AtA)

; 

; Tr: trace of (AtA)^-1.

;

; PROCEDURE:

; Calculate generalized inverse of A by SVD.

;

; MODIFICATION HISTORY:

;   Written by: L.S. (UniBo) 2010 

;      as part of MODAL_INT_MATRIX and ZONAL_INT_MATRIX

;   1) E. D. (INAF-OABo), April 2011: defined as stand-alone routine.

;-

function svd_inv_matrix, a, TOL = tol, DOUBLE = double, sv, tr, DIAG = diag, _EXTRA = extra

on_error, 2

if n_elements(tol) eq 0 then tol = (machar(DOUBLE = double)).eps

; SVD

svdc, transpose(a) ## a, sv, u, v, DOUBLE = double, _EXTRA = extra

; Invert singular values

sv_inv = sv

sv_min = max(abs(sv)) * tol

w = where(abs(sv) ge sv_min, n)  &  if  n ne 0  then  sv_inv[w] = 1 / sv[w]

w = where(abs(sv) lt sv_min, n)  &  if  n ne 0  then  sv_inv[w] = 0

; Invert matrix

n = n_elements(sv)

id = identity(n, DOUBLE = double)

s = lindgen(n)

id[s,s] = sv_inv

b = transpose(temporary(v)) # temporary(id) # temporary(u)

;b = transpose(v) # id # u

tr = trace(b)

diag = diag_matrix(b)

b = temporary(b) ## transpose(a)

return, b

end