00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #if defined (__GNUG__) && defined (USE_PRAGMA_INTERFACE_IMPLEMENTATION)
00024 #pragma implementation
00025 #endif
00026
00027 #ifdef HAVE_CONFIG_H
00028 #include <config.h>
00029 #endif
00030
00031 #include <iostream>
00032
00033 #include "dbleSVD.h"
00034 #include "f77-fcn.h"
00035
00036 extern "C"
00037 {
00038 F77_RET_T
00039 F77_FUNC (dgesvd, DGESVD) (F77_CONST_CHAR_ARG_DECL,
00040 F77_CONST_CHAR_ARG_DECL,
00041 const int&, const int&, double*,
00042 const int&, double*, double*,
00043 const int&, double*, const int&,
00044 double*, const int&, int&
00045 F77_CHAR_ARG_LEN_DECL
00046 F77_CHAR_ARG_LEN_DECL);
00047 }
00048
00049 Matrix
00050 SVD::left_singular_matrix (void) const
00051 {
00052 if (type_computed == SVD::sigma_only)
00053 {
00054 (*current_liboctave_error_handler)
00055 ("ComplexSVD: U not computed because type == SVD::sigma_only");
00056 return Matrix ();
00057 }
00058 else
00059 return left_sm;
00060 }
00061
00062 Matrix
00063 SVD::right_singular_matrix (void) const
00064 {
00065 if (type_computed == SVD::sigma_only)
00066 {
00067 (*current_liboctave_error_handler)
00068 ("ComplexSVD: V not computed because type == SVD::sigma_only");
00069 return Matrix ();
00070 }
00071 else
00072 return right_sm;
00073 }
00074
00075 int
00076 SVD::init (const Matrix& a, SVD::type svd_type)
00077 {
00078 int info;
00079
00080 int m = a.rows ();
00081 int n = a.cols ();
00082
00083 Matrix atmp = a;
00084 double *tmp_data = atmp.fortran_vec ();
00085
00086 int min_mn = m < n ? m : n;
00087
00088 char jobu = 'A';
00089 char jobv = 'A';
00090
00091 int ncol_u = m;
00092 int nrow_vt = n;
00093 int nrow_s = m;
00094 int ncol_s = n;
00095
00096 switch (svd_type)
00097 {
00098 case SVD::economy:
00099 jobu = jobv = 'S';
00100 ncol_u = nrow_vt = nrow_s = ncol_s = min_mn;
00101 break;
00102
00103 case SVD::sigma_only:
00104
00105
00106
00107
00108
00109
00110
00111
00112
00113 jobu = 'N';
00114 jobv = 'N';
00115 ncol_u = nrow_vt = 1;
00116 break;
00117
00118 default:
00119 break;
00120 }
00121
00122 type_computed = svd_type;
00123
00124 if (! (jobu == 'N' || jobu == 'O'))
00125 left_sm.resize (m, ncol_u);
00126
00127 double *u = left_sm.fortran_vec ();
00128
00129 sigma.resize (nrow_s, ncol_s);
00130 double *s_vec = sigma.fortran_vec ();
00131
00132 if (! (jobv == 'N' || jobv == 'O'))
00133 right_sm.resize (nrow_vt, n);
00134
00135 double *vt = right_sm.fortran_vec ();
00136
00137
00138
00139 int lwork = -1;
00140
00141 Array<double> work (1);
00142
00143 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00144 F77_CONST_CHAR_ARG2 (&jobv, 1),
00145 m, n, tmp_data, m, s_vec, u, m, vt,
00146 nrow_vt, work.fortran_vec (), lwork, info
00147 F77_CHAR_ARG_LEN (1)
00148 F77_CHAR_ARG_LEN (1)));
00149
00150 if (f77_exception_encountered)
00151 (*current_liboctave_error_handler) ("unrecoverable error in dgesvd");
00152 else
00153 {
00154 lwork = static_cast<int> (work(0));
00155 work.resize (lwork);
00156
00157 F77_XFCN (dgesvd, DGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00158 F77_CONST_CHAR_ARG2 (&jobv, 1),
00159 m, n, tmp_data, m, s_vec, u, m, vt,
00160 nrow_vt, work.fortran_vec (), lwork, info
00161 F77_CHAR_ARG_LEN (1)
00162 F77_CHAR_ARG_LEN (1)));
00163
00164 if (f77_exception_encountered)
00165 (*current_liboctave_error_handler) ("unrecoverable error in dgesvd");
00166 else
00167 {
00168 if (! (jobv == 'N' || jobv == 'O'))
00169 right_sm = right_sm.transpose ();
00170 }
00171 }
00172
00173 return info;
00174 }
00175
00176 std::ostream&
00177 operator << (std::ostream& os, const SVD& a)
00178 {
00179 os << a.left_singular_matrix () << "\n";
00180 os << a.singular_values () << "\n";
00181 os << a.right_singular_matrix () << "\n";
00182
00183 return os;
00184 }
00185
00186
00187
00188
00189
00190