diff --git a/svd_inv_matrix.pro b/svd_inv_matrix.pro
new file mode 100644
index 0000000000000000000000000000000000000000..f68713e9183d0a2d3f727383f7c38cc74209a746
--- /dev/null
+++ b/svd_inv_matrix.pro
@@ -0,0 +1,70 @@
+; $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