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 "CmplxSVD.h"
00032 #include "f77-fcn.h"
00033 #include "lo-error.h"
00034
00035 extern "C"
00036 {
00037 F77_RET_T
00038 F77_FUNC (zgesvd, ZGESVD) (F77_CONST_CHAR_ARG_DECL,
00039 F77_CONST_CHAR_ARG_DECL,
00040 const int&, const int&, Complex*,
00041 const int&, double*, Complex*, const int&,
00042 Complex*, const int&, Complex*, const int&,
00043 double*, int&
00044 F77_CHAR_ARG_LEN_DECL
00045 F77_CHAR_ARG_LEN_DECL);
00046 }
00047
00048 ComplexMatrix
00049 ComplexSVD::left_singular_matrix (void) const
00050 {
00051 if (type_computed == SVD::sigma_only)
00052 {
00053 (*current_liboctave_error_handler)
00054 ("ComplexSVD: U not computed because type == SVD::sigma_only");
00055 return ComplexMatrix ();
00056 }
00057 else
00058 return left_sm;
00059 }
00060
00061 ComplexMatrix
00062 ComplexSVD::right_singular_matrix (void) const
00063 {
00064 if (type_computed == SVD::sigma_only)
00065 {
00066 (*current_liboctave_error_handler)
00067 ("ComplexSVD: V not computed because type == SVD::sigma_only");
00068 return ComplexMatrix ();
00069 }
00070 else
00071 return right_sm;
00072 }
00073
00074 int
00075 ComplexSVD::init (const ComplexMatrix& a, SVD::type svd_type)
00076 {
00077 int info;
00078
00079 int m = a.rows ();
00080 int n = a.cols ();
00081
00082 ComplexMatrix atmp = a;
00083 Complex *tmp_data = atmp.fortran_vec ();
00084
00085 int min_mn = m < n ? m : n;
00086 int max_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 Complex *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 Complex *vt = right_sm.fortran_vec ();
00136
00137 int lrwork = 5*max_mn;
00138
00139 Array<double> rwork (lrwork);
00140
00141
00142
00143 int lwork = -1;
00144
00145 Array<Complex> work (1);
00146
00147 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00148 F77_CONST_CHAR_ARG2 (&jobv, 1),
00149 m, n, tmp_data, m, s_vec, u, m, vt,
00150 nrow_vt, work.fortran_vec (), lwork,
00151 rwork.fortran_vec (), info
00152 F77_CHAR_ARG_LEN (1)
00153 F77_CHAR_ARG_LEN (1)));
00154
00155 if (f77_exception_encountered)
00156 (*current_liboctave_error_handler) ("unrecoverable error in zgesvd");
00157 else
00158 {
00159 lwork = static_cast<int> (work(0).real ());
00160 work.resize (lwork);
00161
00162 F77_XFCN (zgesvd, ZGESVD, (F77_CONST_CHAR_ARG2 (&jobu, 1),
00163 F77_CONST_CHAR_ARG2 (&jobv, 1),
00164 m, n, tmp_data, m, s_vec, u, m, vt,
00165 nrow_vt, work.fortran_vec (), lwork,
00166 rwork.fortran_vec (), info
00167 F77_CHAR_ARG_LEN (1)
00168 F77_CHAR_ARG_LEN (1)));
00169
00170 if (f77_exception_encountered)
00171 (*current_liboctave_error_handler) ("unrecoverable error in zgesvd");
00172 else
00173 {
00174 if (! (jobv == 'N' || jobv == 'O'))
00175 right_sm = right_sm.hermitian ();
00176 }
00177 }
00178
00179 return info;
00180 }
00181
00182
00183
00184
00185
00186