From 014e9147aac9c8ac664967c80ac2b4f624ce9eb7 Mon Sep 17 00:00:00 2001
From: Guido Kanschat <dr.guido.kanschat@gmail.com>
Date: Fri, 5 Nov 2010 15:30:06 +0000
Subject: [PATCH] start some work on SVD

git-svn-id: https://svn.dealii.org/trunk@22611 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/configure.in                           |  2 +-
 .../include/deal.II/lac/lapack_full_matrix.h   | 18 +++++++++++++++---
 deal.II/include/deal.II/lac/lapack_support.h   |  2 ++
 .../include/deal.II/lac/lapack_templates.h.in  |  9 +++++++++
 4 files changed, 27 insertions(+), 4 deletions(-)

diff --git a/deal.II/configure.in b/deal.II/configure.in
index 7dab4cc12b..f1a5b5c1ab 100644
--- a/deal.II/configure.in
+++ b/deal.II/configure.in
@@ -621,7 +621,7 @@ dnl exception.
 AC_CHECK_FUNCS([daxpy_ saxpy_ dgemv_ sgemv_ dgeev_ sgeev_ dgeevx_ sgeevx_])
 AC_CHECK_FUNCS([dgesvd_ sgesvd_ dgemm_ sgemm_ dgetrf_ sgetrf_ dgetri_ sgetri_])
 AC_CHECK_FUNCS([dgeqrf_ sgeqrf_ dormqr_ sormqr_ dorgqr_ sorgqr_ dtrtrs_ strtrs_])
-AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_])
+AC_CHECK_FUNCS([dgetrs_ sgetrs_ dstev_ sstev_ dsygv_ ssygv_ dgelsd_ sgelsd_])
 
 dnl -------------------------------------------------------------
 dnl       set include paths of several libraries
diff --git a/deal.II/include/deal.II/lac/lapack_full_matrix.h b/deal.II/include/deal.II/lac/lapack_full_matrix.h
index 1d88428a00..52506732b0 100644
--- a/deal.II/include/deal.II/lac/lapack_full_matrix.h
+++ b/deal.II/include/deal.II/lac/lapack_full_matrix.h
@@ -20,6 +20,7 @@
 #include <lac/lapack_support.h>
 #include <lac/vector_memory.h>
 
+#include <base/std_cxx1x/shared_ptr.h>
 #include <vector>
 #include <complex>
 
@@ -436,9 +437,9 @@ class LAPACKFullMatrix : public TransposeTable<number>
     std::vector<number> inv_work;
 
 				     /**
-				      * Real parts of
-				      * eigenvalues. Filled by
-				      * compute_eigenvalues.
+				      * Real parts of eigenvalues or
+				      * the singular values. Filled by
+				      * compute_eigenvalues() or compute_svd().
 				      */
     std::vector<number> wr;
 
@@ -460,6 +461,17 @@ class LAPACKFullMatrix : public TransposeTable<number>
 				      * can be stored.
 				      */
     std::vector<number> vr;
+    
+    /**
+     * The matrix <i>U</i> in the singular value decomposition
+     * <i>USV<sup>T</sup></i>.
+     */
+    boost::shared_ptr<LAPACKFullMatrix<number> > svd_u;
+    /**
+     * The matrix <i>V<sup>T</sup></i> in the singular value decomposition
+     * <i>USV<sup>T</sup></i>.
+     */
+    boost::shared_ptr<LAPACKFullMatrix<number> > svd_vt;
 };
 
 
diff --git a/deal.II/include/deal.II/lac/lapack_support.h b/deal.II/include/deal.II/lac/lapack_support.h
index b23eaa14f5..6aa4b9ee05 100644
--- a/deal.II/include/deal.II/lac/lapack_support.h
+++ b/deal.II/include/deal.II/lac/lapack_support.h
@@ -40,6 +40,8 @@ namespace LAPACKSupport
 	lu,
 					 /// Eigenvalue vector is filled
 	eigenvalues,
+	                                 /// Matrix contains singular value decomposition,
+	svd,
 					 /// Contents is something useless.
 	unusable = 0x8000
   };
diff --git a/deal.II/include/deal.II/lac/lapack_templates.h.in b/deal.II/include/deal.II/lac/lapack_templates.h.in
index 8ddee8562c..631796654d 100644
--- a/deal.II/include/deal.II/lac/lapack_templates.h.in
+++ b/deal.II/include/deal.II/lac/lapack_templates.h.in
@@ -77,6 +77,7 @@ void dsygv_ (const int* itype, const char* jobz, const char* uplo,
              const int* n, double* A, const int* lda, double* B,
              const int* ldb, double* w, double* work,
              const int* lwork, int* info);
+
 // Compute singular value decomposition
 void dgesvd_ (int* jobu, int* jobvt,
 	      const int* n, const int* m, double* A, const int* lda,
@@ -86,6 +87,14 @@ void dgesvd_ (int* jobu, int* jobvt,
 	      double* work, const int* lwork,
 	      int* info);
 
+// Solve a least squares problem using SVD
+void dgelsd_ (const int* m, const int* n, const int* nrhs,
+	      const double* A, const int* lda,
+	      double* B, const int* ldb,
+	      double* s, const double* rcond,
+	      int* rank,
+	      double* work, const int* lwork, int* iwork,
+	      int* info);
 
 // Symmetric tridiagonal matrix
 void dstev_ (const char* jobz, const int* n,
-- 
2.39.5