diff --git a/dspl/src/blas.h b/dspl/src/blas.h
index b0b67d6..8175075 100644
--- a/dspl/src/blas.h
+++ b/dspl/src/blas.h
@@ -362,5 +362,10 @@ int FORTRAN_FUNC(xher2m)(const char*, const char*, const char*, const int*, cons
 
 void FORTRAN_FUNC(zgees)(const char*, const char*, void*, int*, complex_t*, int*, int*, complex_t*, complex_t*, int*, complex_t*, int*, double*, int*, int*);
 
+void FORTRAN_FUNC(dgesdd)(const char*, int*, int*, double*, int*, double*, 
+                          double*, int*, double*, int*, double*, 
+                          int*, int*, int*);
+
+
 #endif
 
diff --git a/dspl/src/matrix.c b/dspl/src/matrix.c
index 2cd2023..83b2d78 100644
--- a/dspl/src/matrix.c
+++ b/dspl/src/matrix.c
@@ -199,6 +199,8 @@ int DSPL_API matrix_eye_cmplx(complex_t* a, int n, int m)
 
 
 
+
+
 #ifdef DOXYGEN_ENGLISH
 
 #endif
@@ -226,6 +228,57 @@ int DSPL_API matrix_mul(double* a, int na, int ma,
 
 
 
+
+#ifdef DOXYGEN_ENGLISH
+
+#endif
+#ifdef DOXYGEN_RUSSIAN
+
+#endif
+int DSPL_API matrix_svd(double* a, int n, int m, 
+                        double* u, double* s, double* vt, int* info)
+{
+    char jobz = 'A';
+    double* work = NULL;
+    int* iwork = NULL;
+    int lwork, mn, mx, err;
+    int pi;
+    
+    if(!a || !u || !s || !vt)
+        return ERROR_PTR;
+    if(n < 1 || m < 1)
+        return ERROR_SIZE;
+      
+    
+    if(n > m)
+    {
+        mn = m;
+        mx = n;
+    }
+    else
+    { 
+        mn = n; 
+        mx = m;
+    }
+    
+    err = RES_OK;
+    
+    lwork = 4 * mn * mn + 6 * mn + mx;
+    work = (double*) malloc(lwork*sizeof(double));
+    iwork = (int*) malloc(8*mn*sizeof(int));
+    dgesdd_(&jobz, &n, &m, a, &n, s, u, &n, vt, &m, work, &lwork, iwork, &pi);
+    
+    if(info)
+        *info = pi;
+    if(pi)
+        err = ERROR_LAPACK;
+    if(work)
+        free(work);
+    if(iwork)
+        free(iwork);
+    return err;
+}
+
 #ifdef DOXYGEN_ENGLISH
 
 #endif
diff --git a/examples/src/matrix_svd_test.c b/examples/src/matrix_svd_test.c
new file mode 100644
index 0000000..88c81db
--- /dev/null
+++ b/examples/src/matrix_svd_test.c
@@ -0,0 +1,65 @@
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include "dspl.h"
+
+#define N  4
+#define M  3
+
+
+int main(int argc, char* argv[])
+{
+    void* handle;           /* DSPL handle         */
+    handle = dspl_load();   /* Load DSPL function  */
+
+    /* Matrix 
+    A = [ 1   2   3;
+          3   4   5;
+          2   0   1;
+          2  -1   0;];
+    in array a by columns
+    */
+    double a[N*M] = { 1, 3, 2, 2,   2, 4, 0, -1,   3, 5, 1, 0};
+    double u[N*N]; /* left orthogonal matrix U              */
+    double s[M];   /* Singular values diagonal              */
+    double vt[M*M];/* transposed left orthogonal matrix V^T */
+
+    double ur[N*M] = {0}; /* matrix UR = U*S */
+    double ar[N*M];       /* AR = UR * V^T   */ 
+    
+    int err, info, i, j, mn;
+    
+    /* print input matrix */
+    matrix_print(a, N, M, "A", "%8.2f");
+
+    /* SVD decomposition A = U*S*V^T */
+    /*-----------------------------------------------------*/
+    err = matrix_svd(a, N, M, u, s, vt, &info);
+    if(err != RES_OK)
+        printf("err = %.8x  info = %d\n", err, info);
+    
+    /* Print SVD decomposition */
+    matrix_print(u,  N, N, "U",   "%8.4f");
+    matrix_print(vt, M, M, "V^T", "%8.4f");
+    matrix_print(s,  M, 1, "S",   "%8.6f");
+
+
+    /* SVD reconstruction AR = U*S*V^T */
+    /*-----------------------------------------------------*/
+
+    /* step 1: UR = U*S. 
+       We can process this by columns 
+       because S is diagonal 
+    */
+    mn = (M > N) ? N : M;
+    for(i = 0; i < mn; i++)
+        for(j = 0; j < N; j++)
+            ur[j + i*N] = u[j + i*N] * s[i];
+
+    /* step 2: AR = U * S * V^T */
+    matrix_mul(ur, N, M, vt, M, M, ar);
+    matrix_print(ar, N, M, "AR = U*S*V^T",   "%8.2f");
+
+    dspl_free(handle);      /* free dspl handle  */
+    return 0;
+}
diff --git a/include/dspl.c b/include/dspl.c
index 1d7a862..b7a818e 100644
--- a/include/dspl.c
+++ b/include/dspl.c
@@ -139,6 +139,7 @@ p_matrix_eye_cmplx                      matrix_eye_cmplx              ;
 p_matrix_mul                            matrix_mul                    ;
 p_matrix_print                          matrix_print                  ;
 p_matrix_print_cmplx                    matrix_print_cmplx            ;
+p_matrix_svd                            matrix_svd                    ;
 p_matrix_transpose                      matrix_transpose              ;
 p_matrix_transpose_cmplx                matrix_transpose_cmplx        ;
 p_matrix_transpose_hermite              matrix_transpose_hermite      ;
@@ -351,6 +352,7 @@ void* dspl_load()
     LOAD_FUNC(matrix_mul);
     LOAD_FUNC(matrix_print);
     LOAD_FUNC(matrix_print_cmplx);
+    LOAD_FUNC(matrix_svd);
     LOAD_FUNC(matrix_transpose);
     LOAD_FUNC(matrix_transpose_cmplx);
     LOAD_FUNC(matrix_transpose_hermite);
@@ -440,8 +442,6 @@ void* dspl_load()
 
 
 
-
-
 void dspl_free(void* handle)
 {
     #ifdef WIN_OS
diff --git a/include/dspl.h b/include/dspl.h
index 09c4d97..63b8489 100644
--- a/include/dspl.h
+++ b/include/dspl.h
@@ -1228,6 +1228,14 @@ DECLARE_FUNC(int,        matrix_print_cmplx,          complex_t*       a
                                                 COMMA const char*      name
                                                 COMMA const char*      format);
 /*----------------------------------------------------------------------------*/
+DECLARE_FUNC(int,        matrix_svd,                  double*          a
+                                                COMMA int              n
+                                                COMMA int              m
+                                                COMMA double*          u
+                                                COMMA double*          s
+                                                COMMA double*          vt
+                                                COMMA int*             info);
+/*----------------------------------------------------------------------------*/
 DECLARE_FUNC(int,        matrix_transpose,            double*          a
                                                 COMMA int              n
                                                 COMMA int              m
@@ -1563,49 +1571,6 @@ DECLARE_FUNC(int,        xcorr_cmplx,                 complex_t*        x
 #endif
 
 
-
-
-
-
-#ifdef DOXYGEN_ENGLISH
-/*! ****************************************************************************
-\ingroup SYS_LOADING_GROUP
-\fn void dspl_free(void* handle)
-\brief Cleans up the previously linked DSPL-2.0 dynamic library.
-
-This cross-platform function clears the library `libdspl.dll` in
-Windows system and from the library `libdspl.so` on the Linux system.
-After cleaning the library, all functions will become unavailable.
-
-\param [in] handle
-Handle of the previously linked DSPL-2.0 library. \n
-This pointer can be `NULL`, in this case no action
-are being produced.
-
-\author Bakhurin Sergey. www.dsplib.org
-***************************************************************************** */
-#endif
-#ifdef DOXYGEN_RUSSIAN
-/*! ****************************************************************************
-\ingroup SYS_LOADING_GROUP
-\fn void dspl_free(void* handle)
-\brief Очищает связанную ранее динамическую библиотеку DSPL-2.0.
-
-Данная кроссплатформенная функция производит очистку библиотеки `libdspl.dll` в 
-системе Windows и с библиотеки `libdspl.so` в системе Linux. 
-После очистки библиотеки все функции станут недоступны.
-
-\param[in] handle
-Хэндл прилинкованной ранее библиотеки DSPL-2.0. \n
-Данный указатель может быть `NULL`, в этом случае никакие действия не 
-производятся.\n\n
-
-\author Бахурин Сергей. www.dsplib.org
-**************************************************************************** */
-#endif
-void* dspl_load();
-
-
 #ifdef DOXYGEN_ENGLISH
 /*! ****************************************************************************
 \ingroup SYS_LOADING_GROUP
@@ -1694,7 +1659,7 @@ int main(int argc, char* argv[])
     void* hdspl;           // DSPL хэндл
     hdspl = dspl_load();   // Динамическая линковка
     
-    // Проверяем указатель. Если `NULLL`, то линковка прошла неудачно
+    // Проверяем указатель. Если `NULL`, то линковка прошла неудачно
     if(!hdspl)
     {
         printf("libdspl loading error!\n");
@@ -1714,6 +1679,48 @@ int main(int argc, char* argv[])
 \author Бахурин Сергей. www.dsplib.org
 ***************************************************************************** */
 #endif
+void* dspl_load();
+
+
+
+
+
+#ifdef DOXYGEN_ENGLISH
+/*! ****************************************************************************
+\ingroup SYS_LOADING_GROUP
+\fn void dspl_free(void* handle)
+\brief Cleans up the previously linked DSPL-2.0 dynamic library.
+
+This cross-platform function clears the library `libdspl.dll` in
+Windows system and from the library `libdspl.so` on the Linux system.
+After cleaning the library, all functions will become unavailable.
+
+\param [in] handle
+Handle of the previously linked DSPL-2.0 library. \n
+This pointer can be `NULL`, in this case no action
+are being produced.
+
+\author Bakhurin Sergey. www.dsplib.org
+***************************************************************************** */
+#endif
+#ifdef DOXYGEN_RUSSIAN
+/*! ****************************************************************************
+\ingroup SYS_LOADING_GROUP
+\fn void dspl_free(void* handle)
+\brief Очищает связанную ранее динамическую библиотеку DSPL-2.0.
+
+Данная кроссплатформенная функция производит очистку библиотеки `libdspl.dll` в 
+системе Windows и с библиотеки `libdspl.so` в системе Linux. 
+После очистки библиотеки все функции станут недоступны.
+
+\param[in] handle
+Хэндл прилинкованной ранее библиотеки DSPL-2.0. \n
+Данный указатель может быть `NULL`, в этом случае никакие действия не 
+производятся.\n\n
+
+\author Бахурин Сергей. www.dsplib.org
+**************************************************************************** */
+#endif
 void  dspl_free(void* handle);