00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 template<typename eT>
00023 inline
00024 void
00025 auxlib::inv_noalias(Mat<eT>& out, const Mat<eT>& X)
00026 {
00027 arma_extra_debug_sigprint();
00028
00029 switch(X.n_rows)
00030 {
00031 case 1:
00032 {
00033 out.set_size(1,1);
00034 out[0] = eT(1) / X[0];
00035 };
00036 break;
00037
00038 case 2:
00039 {
00040 out.set_size(2,2);
00041
00042 const eT a = X.at(0,0);
00043 const eT b = X.at(0,1);
00044 const eT c = X.at(1,0);
00045 const eT d = X.at(1,1);
00046
00047 const eT k = eT(1) / (a*d - b*c);
00048
00049 out.at(0,0) = d*k;
00050 out.at(0,1) = -b*k;
00051 out.at(1,0) = -c*k;
00052 out.at(1,1) = a*k;
00053 };
00054 break;
00055
00056 case 3:
00057 {
00058 out.set_size(3,3);
00059
00060 const eT* X_col0 = X.colptr(0);
00061 const eT a11 = X_col0[0];
00062 const eT a21 = X_col0[1];
00063 const eT a31 = X_col0[2];
00064
00065 const eT* X_col1 = X.colptr(1);
00066 const eT a12 = X_col1[0];
00067 const eT a22 = X_col1[1];
00068 const eT a32 = X_col1[2];
00069
00070 const eT* X_col2 = X.colptr(2);
00071 const eT a13 = X_col2[0];
00072 const eT a23 = X_col2[1];
00073 const eT a33 = X_col2[2];
00074
00075 const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00076
00077
00078 eT* out_col0 = out.colptr(0);
00079 out_col0[0] = (a33*a22 - a32*a23) * k;
00080 out_col0[1] = -(a33*a21 - a31*a23) * k;
00081 out_col0[2] = (a32*a21 - a31*a22) * k;
00082
00083 eT* out_col1 = out.colptr(1);
00084 out_col1[0] = -(a33*a12 - a32*a13) * k;
00085 out_col1[1] = (a33*a11 - a31*a13) * k;
00086 out_col1[2] = -(a32*a11 - a31*a12) * k;
00087
00088 eT* out_col2 = out.colptr(2);
00089 out_col2[0] = (a23*a12 - a22*a13) * k;
00090 out_col2[1] = -(a23*a11 - a21*a13) * k;
00091 out_col2[2] = (a22*a11 - a21*a12) * k;
00092 };
00093 break;
00094
00095 case 4:
00096 {
00097 out.set_size(4,4);
00098
00099 out.at(0,0) = X.at(1,2)*X.at(2,3)*X.at(3,1) - X.at(1,3)*X.at(2,2)*X.at(3,1) + X.at(1,3)*X.at(2,1)*X.at(3,2) - X.at(1,1)*X.at(2,3)*X.at(3,2) - X.at(1,2)*X.at(2,1)*X.at(3,3) + X.at(1,1)*X.at(2,2)*X.at(3,3);
00100 out.at(1,0) = X.at(1,3)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,0)*X.at(3,2) + X.at(1,0)*X.at(2,3)*X.at(3,2) + X.at(1,2)*X.at(2,0)*X.at(3,3) - X.at(1,0)*X.at(2,2)*X.at(3,3);
00101 out.at(2,0) = X.at(1,1)*X.at(2,3)*X.at(3,0) - X.at(1,3)*X.at(2,1)*X.at(3,0) + X.at(1,3)*X.at(2,0)*X.at(3,1) - X.at(1,0)*X.at(2,3)*X.at(3,1) - X.at(1,1)*X.at(2,0)*X.at(3,3) + X.at(1,0)*X.at(2,1)*X.at(3,3);
00102 out.at(3,0) = X.at(1,2)*X.at(2,1)*X.at(3,0) - X.at(1,1)*X.at(2,2)*X.at(3,0) - X.at(1,2)*X.at(2,0)*X.at(3,1) + X.at(1,0)*X.at(2,2)*X.at(3,1) + X.at(1,1)*X.at(2,0)*X.at(3,2) - X.at(1,0)*X.at(2,1)*X.at(3,2);
00103
00104 out.at(0,1) = X.at(0,3)*X.at(2,2)*X.at(3,1) - X.at(0,2)*X.at(2,3)*X.at(3,1) - X.at(0,3)*X.at(2,1)*X.at(3,2) + X.at(0,1)*X.at(2,3)*X.at(3,2) + X.at(0,2)*X.at(2,1)*X.at(3,3) - X.at(0,1)*X.at(2,2)*X.at(3,3);
00105 out.at(1,1) = X.at(0,2)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,2)*X.at(3,0) + X.at(0,3)*X.at(2,0)*X.at(3,2) - X.at(0,0)*X.at(2,3)*X.at(3,2) - X.at(0,2)*X.at(2,0)*X.at(3,3) + X.at(0,0)*X.at(2,2)*X.at(3,3);
00106 out.at(2,1) = X.at(0,3)*X.at(2,1)*X.at(3,0) - X.at(0,1)*X.at(2,3)*X.at(3,0) - X.at(0,3)*X.at(2,0)*X.at(3,1) + X.at(0,0)*X.at(2,3)*X.at(3,1) + X.at(0,1)*X.at(2,0)*X.at(3,3) - X.at(0,0)*X.at(2,1)*X.at(3,3);
00107 out.at(3,1) = X.at(0,1)*X.at(2,2)*X.at(3,0) - X.at(0,2)*X.at(2,1)*X.at(3,0) + X.at(0,2)*X.at(2,0)*X.at(3,1) - X.at(0,0)*X.at(2,2)*X.at(3,1) - X.at(0,1)*X.at(2,0)*X.at(3,2) + X.at(0,0)*X.at(2,1)*X.at(3,2);
00108
00109 out.at(0,2) = X.at(0,2)*X.at(1,3)*X.at(3,1) - X.at(0,3)*X.at(1,2)*X.at(3,1) + X.at(0,3)*X.at(1,1)*X.at(3,2) - X.at(0,1)*X.at(1,3)*X.at(3,2) - X.at(0,2)*X.at(1,1)*X.at(3,3) + X.at(0,1)*X.at(1,2)*X.at(3,3);
00110 out.at(1,2) = X.at(0,3)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,0)*X.at(3,2) + X.at(0,0)*X.at(1,3)*X.at(3,2) + X.at(0,2)*X.at(1,0)*X.at(3,3) - X.at(0,0)*X.at(1,2)*X.at(3,3);
00111 out.at(2,2) = X.at(0,1)*X.at(1,3)*X.at(3,0) - X.at(0,3)*X.at(1,1)*X.at(3,0) + X.at(0,3)*X.at(1,0)*X.at(3,1) - X.at(0,0)*X.at(1,3)*X.at(3,1) - X.at(0,1)*X.at(1,0)*X.at(3,3) + X.at(0,0)*X.at(1,1)*X.at(3,3);
00112 out.at(3,2) = X.at(0,2)*X.at(1,1)*X.at(3,0) - X.at(0,1)*X.at(1,2)*X.at(3,0) - X.at(0,2)*X.at(1,0)*X.at(3,1) + X.at(0,0)*X.at(1,2)*X.at(3,1) + X.at(0,1)*X.at(1,0)*X.at(3,2) - X.at(0,0)*X.at(1,1)*X.at(3,2);
00113
00114 out.at(0,3) = X.at(0,3)*X.at(1,2)*X.at(2,1) - X.at(0,2)*X.at(1,3)*X.at(2,1) - X.at(0,3)*X.at(1,1)*X.at(2,2) + X.at(0,1)*X.at(1,3)*X.at(2,2) + X.at(0,2)*X.at(1,1)*X.at(2,3) - X.at(0,1)*X.at(1,2)*X.at(2,3);
00115 out.at(1,3) = X.at(0,2)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,2)*X.at(2,0) + X.at(0,3)*X.at(1,0)*X.at(2,2) - X.at(0,0)*X.at(1,3)*X.at(2,2) - X.at(0,2)*X.at(1,0)*X.at(2,3) + X.at(0,0)*X.at(1,2)*X.at(2,3);
00116 out.at(2,3) = X.at(0,3)*X.at(1,1)*X.at(2,0) - X.at(0,1)*X.at(1,3)*X.at(2,0) - X.at(0,3)*X.at(1,0)*X.at(2,1) + X.at(0,0)*X.at(1,3)*X.at(2,1) + X.at(0,1)*X.at(1,0)*X.at(2,3) - X.at(0,0)*X.at(1,1)*X.at(2,3);
00117 out.at(3,3) = X.at(0,1)*X.at(1,2)*X.at(2,0) - X.at(0,2)*X.at(1,1)*X.at(2,0) + X.at(0,2)*X.at(1,0)*X.at(2,1) - X.at(0,0)*X.at(1,2)*X.at(2,1) - X.at(0,1)*X.at(1,0)*X.at(2,2) + X.at(0,0)*X.at(1,1)*X.at(2,2);
00118
00119 out /= det(X);
00120 };
00121 break;
00122
00123 default:
00124 {
00125 #if defined(ARMA_USE_ATLAS)
00126 {
00127 out = X;
00128 podarray<int> ipiv(out.n_rows);
00129
00130 int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
00131
00132 if(info == 0)
00133 {
00134 info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
00135 }
00136
00137 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00138 }
00139 #elif defined(ARMA_USE_LAPACK)
00140 {
00141 out = X;
00142
00143 int n_rows = out.n_rows;
00144 int n_cols = out.n_cols;
00145 int info = 0;
00146
00147 podarray<int> ipiv(out.n_rows);
00148
00149
00150
00151
00152
00153
00154
00155 int work_len = (std::max)(1, n_rows*84);
00156 podarray<eT> work(work_len);
00157
00158 lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
00159
00160 if(info == 0)
00161 {
00162
00163
00164 int work_len_tmp = -1;
00165 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
00166
00167 if(info == 0)
00168 {
00169 int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0]));
00170
00171
00172 if(work_len < proposed_work_len)
00173 {
00174 work_len = proposed_work_len;
00175 work.set_size(work_len);
00176 }
00177 }
00178
00179 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
00180 }
00181
00182 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00183 }
00184 #else
00185 {
00186 arma_stop("auxlib::inv_noalias(): need ATLAS or LAPACK library");
00187 }
00188 #endif
00189 };
00190 }
00191 }
00192
00193
00194
00195
00196 template<typename eT>
00197 inline
00198 void
00199 auxlib::inv_inplace(Mat<eT>& X)
00200 {
00201 arma_extra_debug_sigprint();
00202
00203
00204
00205
00206
00207
00208
00209 switch(X.n_rows)
00210 {
00211 case 1:
00212 {
00213 X[0] = eT(1) / X[0];
00214 };
00215 break;
00216
00217 case 2:
00218 {
00219 const eT a = X.at(0,0);
00220 const eT b = X.at(0,1);
00221 const eT c = X.at(1,0);
00222 const eT d = X.at(1,1);
00223
00224 const eT k = eT(1) / (a*d - b*c);
00225
00226 X.at(0,0) = d*k;
00227 X.at(0,1) = -b*k;
00228 X.at(1,0) = -c*k;
00229 X.at(1,1) = a*k;
00230 };
00231 break;
00232
00233 case 3:
00234 {
00235 eT* X_col0 = X.colptr(0);
00236 eT* X_col1 = X.colptr(1);
00237 eT* X_col2 = X.colptr(2);
00238
00239 const eT a11 = X_col0[0];
00240 const eT a21 = X_col0[1];
00241 const eT a31 = X_col0[2];
00242
00243 const eT a12 = X_col1[0];
00244 const eT a22 = X_col1[1];
00245 const eT a32 = X_col1[2];
00246
00247 const eT a13 = X_col2[0];
00248 const eT a23 = X_col2[1];
00249 const eT a33 = X_col2[2];
00250
00251 const eT k = eT(1) / ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00252
00253 X_col0[0] = (a33*a22 - a32*a23) * k;
00254 X_col0[1] = -(a33*a21 - a31*a23) * k;
00255 X_col0[2] = (a32*a21 - a31*a22) * k;
00256
00257 X_col1[0] = -(a33*a12 - a32*a13) * k;
00258 X_col1[1] = (a33*a11 - a31*a13) * k;
00259 X_col1[2] = -(a32*a11 - a31*a12) * k;
00260
00261 X_col2[0] = (a23*a12 - a22*a13) * k;
00262 X_col2[1] = -(a23*a11 - a21*a13) * k;
00263 X_col2[2] = (a22*a11 - a21*a12) * k;
00264 };
00265 break;
00266
00267 case 4:
00268 {
00269 const Mat<eT> A(X);
00270
00271 X.at(0,0) = A.at(1,2)*A.at(2,3)*A.at(3,1) - A.at(1,3)*A.at(2,2)*A.at(3,1) + A.at(1,3)*A.at(2,1)*A.at(3,2) - A.at(1,1)*A.at(2,3)*A.at(3,2) - A.at(1,2)*A.at(2,1)*A.at(3,3) + A.at(1,1)*A.at(2,2)*A.at(3,3);
00272 X.at(1,0) = A.at(1,3)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,0)*A.at(3,2) + A.at(1,0)*A.at(2,3)*A.at(3,2) + A.at(1,2)*A.at(2,0)*A.at(3,3) - A.at(1,0)*A.at(2,2)*A.at(3,3);
00273 X.at(2,0) = A.at(1,1)*A.at(2,3)*A.at(3,0) - A.at(1,3)*A.at(2,1)*A.at(3,0) + A.at(1,3)*A.at(2,0)*A.at(3,1) - A.at(1,0)*A.at(2,3)*A.at(3,1) - A.at(1,1)*A.at(2,0)*A.at(3,3) + A.at(1,0)*A.at(2,1)*A.at(3,3);
00274 X.at(3,0) = A.at(1,2)*A.at(2,1)*A.at(3,0) - A.at(1,1)*A.at(2,2)*A.at(3,0) - A.at(1,2)*A.at(2,0)*A.at(3,1) + A.at(1,0)*A.at(2,2)*A.at(3,1) + A.at(1,1)*A.at(2,0)*A.at(3,2) - A.at(1,0)*A.at(2,1)*A.at(3,2);
00275
00276 X.at(0,1) = A.at(0,3)*A.at(2,2)*A.at(3,1) - A.at(0,2)*A.at(2,3)*A.at(3,1) - A.at(0,3)*A.at(2,1)*A.at(3,2) + A.at(0,1)*A.at(2,3)*A.at(3,2) + A.at(0,2)*A.at(2,1)*A.at(3,3) - A.at(0,1)*A.at(2,2)*A.at(3,3);
00277 X.at(1,1) = A.at(0,2)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,2)*A.at(3,0) + A.at(0,3)*A.at(2,0)*A.at(3,2) - A.at(0,0)*A.at(2,3)*A.at(3,2) - A.at(0,2)*A.at(2,0)*A.at(3,3) + A.at(0,0)*A.at(2,2)*A.at(3,3);
00278 X.at(2,1) = A.at(0,3)*A.at(2,1)*A.at(3,0) - A.at(0,1)*A.at(2,3)*A.at(3,0) - A.at(0,3)*A.at(2,0)*A.at(3,1) + A.at(0,0)*A.at(2,3)*A.at(3,1) + A.at(0,1)*A.at(2,0)*A.at(3,3) - A.at(0,0)*A.at(2,1)*A.at(3,3);
00279 X.at(3,1) = A.at(0,1)*A.at(2,2)*A.at(3,0) - A.at(0,2)*A.at(2,1)*A.at(3,0) + A.at(0,2)*A.at(2,0)*A.at(3,1) - A.at(0,0)*A.at(2,2)*A.at(3,1) - A.at(0,1)*A.at(2,0)*A.at(3,2) + A.at(0,0)*A.at(2,1)*A.at(3,2);
00280
00281 X.at(0,2) = A.at(0,2)*A.at(1,3)*A.at(3,1) - A.at(0,3)*A.at(1,2)*A.at(3,1) + A.at(0,3)*A.at(1,1)*A.at(3,2) - A.at(0,1)*A.at(1,3)*A.at(3,2) - A.at(0,2)*A.at(1,1)*A.at(3,3) + A.at(0,1)*A.at(1,2)*A.at(3,3);
00282 X.at(1,2) = A.at(0,3)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,0)*A.at(3,2) + A.at(0,0)*A.at(1,3)*A.at(3,2) + A.at(0,2)*A.at(1,0)*A.at(3,3) - A.at(0,0)*A.at(1,2)*A.at(3,3);
00283 X.at(2,2) = A.at(0,1)*A.at(1,3)*A.at(3,0) - A.at(0,3)*A.at(1,1)*A.at(3,0) + A.at(0,3)*A.at(1,0)*A.at(3,1) - A.at(0,0)*A.at(1,3)*A.at(3,1) - A.at(0,1)*A.at(1,0)*A.at(3,3) + A.at(0,0)*A.at(1,1)*A.at(3,3);
00284 X.at(3,2) = A.at(0,2)*A.at(1,1)*A.at(3,0) - A.at(0,1)*A.at(1,2)*A.at(3,0) - A.at(0,2)*A.at(1,0)*A.at(3,1) + A.at(0,0)*A.at(1,2)*A.at(3,1) + A.at(0,1)*A.at(1,0)*A.at(3,2) - A.at(0,0)*A.at(1,1)*A.at(3,2);
00285
00286 X.at(0,3) = A.at(0,3)*A.at(1,2)*A.at(2,1) - A.at(0,2)*A.at(1,3)*A.at(2,1) - A.at(0,3)*A.at(1,1)*A.at(2,2) + A.at(0,1)*A.at(1,3)*A.at(2,2) + A.at(0,2)*A.at(1,1)*A.at(2,3) - A.at(0,1)*A.at(1,2)*A.at(2,3);
00287 X.at(1,3) = A.at(0,2)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,2)*A.at(2,0) + A.at(0,3)*A.at(1,0)*A.at(2,2) - A.at(0,0)*A.at(1,3)*A.at(2,2) - A.at(0,2)*A.at(1,0)*A.at(2,3) + A.at(0,0)*A.at(1,2)*A.at(2,3);
00288 X.at(2,3) = A.at(0,3)*A.at(1,1)*A.at(2,0) - A.at(0,1)*A.at(1,3)*A.at(2,0) - A.at(0,3)*A.at(1,0)*A.at(2,1) + A.at(0,0)*A.at(1,3)*A.at(2,1) + A.at(0,1)*A.at(1,0)*A.at(2,3) - A.at(0,0)*A.at(1,1)*A.at(2,3);
00289 X.at(3,3) = A.at(0,1)*A.at(1,2)*A.at(2,0) - A.at(0,2)*A.at(1,1)*A.at(2,0) + A.at(0,2)*A.at(1,0)*A.at(2,1) - A.at(0,0)*A.at(1,2)*A.at(2,1) - A.at(0,1)*A.at(1,0)*A.at(2,2) + A.at(0,0)*A.at(1,1)*A.at(2,2);
00290
00291 X /= det(A);
00292 };
00293 break;
00294
00295 default:
00296 {
00297 #if defined(ARMA_USE_ATLAS)
00298 {
00299 Mat<eT>& out = X;
00300 podarray<int> ipiv(out.n_rows);
00301
00302 int info = atlas::clapack_getrf(atlas::CblasColMajor, out.n_rows, out.n_cols, out.memptr(), out.n_rows, ipiv.memptr());
00303
00304 if(info == 0)
00305 {
00306 info = atlas::clapack_getri(atlas::CblasColMajor, out.n_rows, out.memptr(), out.n_rows, ipiv.memptr());
00307 }
00308
00309 arma_check( (info > 0), "auxlib::inv_inplace(): matrix appears to be singular" );
00310 }
00311 #elif defined(ARMA_USE_LAPACK)
00312 {
00313 Mat<eT>& out = X;
00314
00315 int n_rows = out.n_rows;
00316 int n_cols = out.n_cols;
00317 int info = 0;
00318
00319 podarray<int> ipiv(out.n_rows);
00320
00321
00322
00323
00324
00325
00326
00327 int work_len = (std::max)(1, n_rows*84);
00328 podarray<eT> work(work_len);
00329
00330 lapack::getrf_(&n_rows, &n_cols, out.memptr(), &n_rows, ipiv.memptr(), &info);
00331
00332 if(info == 0)
00333 {
00334
00335
00336 int work_len_tmp = -1;
00337 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len_tmp, &info);
00338
00339 if(info == 0)
00340 {
00341 int proposed_work_len = static_cast<int>(auxlib::tmp_real(work[0]));
00342
00343
00344 if(work_len < proposed_work_len)
00345 {
00346 work_len = proposed_work_len;
00347 work.set_size(work_len);
00348 }
00349 }
00350
00351 lapack::getri_(&n_rows, out.memptr(), &n_rows, ipiv.memptr(), work.memptr(), &work_len, &info);
00352 }
00353
00354 arma_check( (info > 0), "auxlib::inv_noalias(): matrix appears to be singular" );
00355 }
00356 #else
00357 {
00358 arma_stop("auxlib::inv_inplace(): need ATLAS or LAPACK");
00359 }
00360 #endif
00361 }
00362
00363 }
00364
00365 }
00366
00367
00368
00369 template<typename eT>
00370 inline
00371 eT
00372 auxlib::det(const Mat<eT>& X)
00373 {
00374 arma_extra_debug_sigprint();
00375
00376 switch(X.n_rows)
00377 {
00378 case 0:
00379 return 0.0;
00380
00381 case 1:
00382 return X[0];
00383
00384 case 2:
00385 return (X.at(0,0)*X.at(1,1) - X.at(0,1)*X.at(1,0));
00386
00387 case 3:
00388 {
00389 const eT* a_col0 = X.colptr(0);
00390 const eT a11 = a_col0[0];
00391 const eT a21 = a_col0[1];
00392 const eT a31 = a_col0[2];
00393
00394 const eT* a_col1 = X.colptr(1);
00395 const eT a12 = a_col1[0];
00396 const eT a22 = a_col1[1];
00397 const eT a32 = a_col1[2];
00398
00399 const eT* a_col2 = X.colptr(2);
00400 const eT a13 = a_col2[0];
00401 const eT a23 = a_col2[1];
00402 const eT a33 = a_col2[2];
00403
00404 return ( a11*(a33*a22 - a32*a23) - a21*(a33*a12-a32*a13) + a31*(a23*a12 - a22*a13) );
00405
00406
00407
00408
00409
00410
00411
00412
00413 }
00414
00415 case 4:
00416 {
00417 const eT val = \
00418 X.at(0,3) * X.at(1,2) * X.at(2,1) * X.at(3,0) \
00419 - X.at(0,2) * X.at(1,3) * X.at(2,1) * X.at(3,0) \
00420 - X.at(0,3) * X.at(1,1) * X.at(2,2) * X.at(3,0) \
00421 + X.at(0,1) * X.at(1,3) * X.at(2,2) * X.at(3,0) \
00422 + X.at(0,2) * X.at(1,1) * X.at(2,3) * X.at(3,0) \
00423 - X.at(0,1) * X.at(1,2) * X.at(2,3) * X.at(3,0) \
00424 - X.at(0,3) * X.at(1,2) * X.at(2,0) * X.at(3,1) \
00425 + X.at(0,2) * X.at(1,3) * X.at(2,0) * X.at(3,1) \
00426 + X.at(0,3) * X.at(1,0) * X.at(2,2) * X.at(3,1) \
00427 - X.at(0,0) * X.at(1,3) * X.at(2,2) * X.at(3,1) \
00428 - X.at(0,2) * X.at(1,0) * X.at(2,3) * X.at(3,1) \
00429 + X.at(0,0) * X.at(1,2) * X.at(2,3) * X.at(3,1) \
00430 + X.at(0,3) * X.at(1,1) * X.at(2,0) * X.at(3,2) \
00431 - X.at(0,1) * X.at(1,3) * X.at(2,0) * X.at(3,2) \
00432 - X.at(0,3) * X.at(1,0) * X.at(2,1) * X.at(3,2) \
00433 + X.at(0,0) * X.at(1,3) * X.at(2,1) * X.at(3,2) \
00434 + X.at(0,1) * X.at(1,0) * X.at(2,3) * X.at(3,2) \
00435 - X.at(0,0) * X.at(1,1) * X.at(2,3) * X.at(3,2) \
00436 - X.at(0,2) * X.at(1,1) * X.at(2,0) * X.at(3,3) \
00437 + X.at(0,1) * X.at(1,2) * X.at(2,0) * X.at(3,3) \
00438 + X.at(0,2) * X.at(1,0) * X.at(2,1) * X.at(3,3) \
00439 - X.at(0,0) * X.at(1,2) * X.at(2,1) * X.at(3,3) \
00440 - X.at(0,1) * X.at(1,0) * X.at(2,2) * X.at(3,3) \
00441 + X.at(0,0) * X.at(1,1) * X.at(2,2) * X.at(3,3) \
00442 ;
00443
00444 return val;
00445 }
00446
00447 default:
00448 {
00449 #if defined(ARMA_USE_ATLAS)
00450 {
00451 Mat<eT> tmp = X;
00452 podarray<int> ipiv(tmp.n_rows);
00453
00454 atlas::clapack_getrf(atlas::CblasColMajor, tmp.n_rows, tmp.n_cols, tmp.memptr(), tmp.n_rows, ipiv.memptr());
00455
00456
00457 eT val = tmp.at(0,0);
00458 for(u32 i=1; i < tmp.n_rows; ++i)
00459 {
00460 val *= tmp.at(i,i);
00461 }
00462
00463 eT sign = eT(1);
00464 for(u32 i=0; i < tmp.n_rows; ++i)
00465 {
00466 if(s32(i) != ipiv.mem[i])
00467 {
00468 sign *= eT(-1);
00469 }
00470 }
00471
00472 return sign*val;
00473 }
00474 #elif defined(ARMA_USE_LAPACK)
00475 {
00476 Mat<eT> tmp = X;
00477 podarray<int> ipiv(tmp.n_rows);
00478
00479 int info = 0;
00480 int n_rows = tmp.n_rows;
00481 int n_cols = tmp.n_cols;
00482
00483 lapack::getrf_(&n_rows, &n_cols, tmp.memptr(), &n_rows, ipiv.memptr(), &info);
00484
00485
00486 eT val = tmp.at(0,0);
00487 for(u32 i=1; i < tmp.n_rows; ++i)
00488 {
00489 val *= tmp.at(i,i);
00490 }
00491
00492 eT sign = eT(1);
00493 for(u32 i=0; i < tmp.n_rows; ++i)
00494 {
00495 if(s32(i) != ipiv.mem[i])
00496 {
00497 sign *= eT(-1);
00498 }
00499 }
00500
00501 return sign*val;
00502 }
00503 #else
00504 {
00505 arma_stop("auxlib::det(): need ATLAS or LAPACK library");
00506 return eT(0);
00507 }
00508 #endif
00509 }
00510 }
00511 }
00512
00513
00514
00515
00516 template<typename eT>
00517 inline
00518 void
00519 auxlib::lu(Mat<eT>& L, Mat<eT>& U, podarray<int>& ipiv, const Mat<eT>& X_orig)
00520 {
00521 arma_extra_debug_sigprint();
00522
00523 arma_debug_check( (&L == &U), "auxlib::lu(): L and U are the same object");
00524
00525 unwrap_check< Mat<eT> > tmp1(X_orig, U);
00526 const Mat<eT>& X_tmp = tmp1.M;
00527
00528 unwrap_check< Mat<eT> > tmp2(X_tmp, L);
00529 const Mat<eT>& X = tmp2.M;
00530
00531 U = X;
00532
00533 #if defined(ARMA_USE_ATLAS) || defined(ARMA_USE_LAPACK)
00534 {
00535
00536 #if defined(ARMA_USE_ATLAS)
00537 {
00538 ipiv.set_size(U.n_rows);
00539
00540
00541 atlas::clapack_getrf(atlas::CblasColMajor, U.n_rows, U.n_cols, U.memptr(), U.n_rows, ipiv.memptr());
00542 }
00543 #elif defined(ARMA_USE_LAPACK)
00544 {
00545 ipiv.set_size(U.n_rows);
00546 int info = 0;
00547
00548 int n_rows = U.n_rows;
00549 int n_cols = U.n_cols;
00550
00551 lapack::getrf_(&n_rows, &n_cols, U.memptr(), &n_rows, ipiv.memptr(), &info);
00552 }
00553 #endif
00554
00555
00556 L.set_size(U.n_rows, U.n_rows);
00557
00558 for(u32 col=0; col<L.n_cols; ++col)
00559 {
00560
00561 for(u32 row=0; row<col; ++row)
00562 {
00563 L.at(row,col) = eT(0);
00564 }
00565
00566 L.at(col,col) = eT(1);
00567
00568 for(u32 row=col+1; row<L.n_rows; ++row)
00569 {
00570 L.at(row,col) = U.at(row,col);
00571 U.at(row,col) = eT(0);
00572 }
00573
00574 }
00575 }
00576 #else
00577 {
00578 arma_stop("auxlib::lu(): need ATLAS or LAPACK library");
00579 }
00580 #endif
00581
00582 }
00583
00584
00585
00586 template<typename eT>
00587 inline
00588 void
00589 auxlib::lu(Mat<eT>& L, Mat<eT>& U, Mat<eT>& P, const Mat<eT>& X)
00590 {
00591 arma_extra_debug_sigprint();
00592
00593 arma_debug_check( ( (&P == &L) || (&P == &U) ), "auxlib::lu(): P is the same object as L or U" );
00594
00595 podarray<int> ipiv;
00596 auxlib::lu(L, U, ipiv, X);
00597
00598 const u32 n = ipiv.n_elem;
00599
00600 P.set_size(n,n);
00601 Mat<eT> ident = eye(n,n);
00602
00603 for(u32 i=n; i>0; --i)
00604 {
00605 const u32 j = i-1;
00606 const u32 k = ipiv[j]-1;
00607
00608 ident.swap_rows(j,k);
00609
00610 if(i == n)
00611 P = ident;
00612 else
00613 P *= ident;
00614
00615 ident.swap_rows(j,k);
00616 }
00617
00618 }
00619
00620
00621
00622 template<typename eT>
00623 inline
00624 void
00625 auxlib::lu(Mat<eT>& L, Mat<eT>& U, const Mat<eT>& X)
00626 {
00627 arma_extra_debug_sigprint();
00628
00629 podarray<int> ipiv;
00630 auxlib::lu(L, U, ipiv, X);
00631 }
00632
00633
00634
00635
00636 template<typename eT>
00637 inline
00638 void
00639 auxlib::eig_sym(Col<eT>& eigval, const Mat<eT>& A_orig)
00640 {
00641 arma_extra_debug_sigprint();
00642
00643 #if defined(ARMA_USE_LAPACK)
00644 {
00645 const unwrap_check<Mat<eT> > tmp(A_orig, eigval);
00646 const Mat<eT>& A = tmp.M;
00647
00648 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square");
00649
00650
00651
00652
00653 char jobz = 'N';
00654 char uplo = 'U';
00655
00656 int n_rows = A.n_rows;
00657 int lwork = (std::max)(1,3*n_rows-1);
00658
00659 eigval.set_size(n_rows);
00660 podarray<eT> work(lwork);
00661
00662 Mat<eT> A_copy = A;
00663 int info;
00664
00665 arma_extra_debug_print("lapack::syev_()");
00666 lapack::syev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
00667 }
00668 #else
00669 {
00670 arma_stop("auxlib::eig_sym(): need LAPACK library");
00671 }
00672 #endif
00673 }
00674
00675
00676
00677
00678 template<typename T>
00679 inline
00680 void
00681 auxlib::eig_sym(Col<T>& eigval, const Mat< std::complex<T> >& A)
00682 {
00683 arma_extra_debug_sigprint();
00684
00685 typedef typename std::complex<T> eT;
00686
00687 #if defined(ARMA_USE_LAPACK)
00688 {
00689 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian");
00690
00691 char jobz = 'N';
00692 char uplo = 'U';
00693
00694 int n_rows = A.n_rows;
00695 int lda = A.n_rows;
00696 int lwork = (std::max)(1, 2*n_rows - 1);
00697
00698 eigval.set_size(n_rows);
00699
00700 podarray<eT> work(lwork);
00701 podarray<T> rwork( (std::max)(1, 3*n_rows - 2) );
00702
00703 Mat<eT> A_copy = A;
00704 int info;
00705
00706 arma_extra_debug_print("lapack::heev_()");
00707 lapack::heev_(&jobz, &uplo, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
00708 }
00709 #else
00710 {
00711 arma_stop("auxlib::eig_sym(): need LAPACK library");
00712 }
00713 #endif
00714 }
00715
00716
00717
00718
00719 template<typename eT>
00720 inline
00721 void
00722 auxlib::eig_sym(Col<eT>& eigval, Mat<eT>& eigvec, const Mat<eT>& A_orig)
00723 {
00724 arma_extra_debug_sigprint();
00725
00726
00727
00728 #if defined(ARMA_USE_LAPACK)
00729 {
00730 const unwrap_check< Mat<eT> > tmp1(A_orig, eigval);
00731 const Mat<eT>& A_tmp = tmp1.M;
00732
00733 const unwrap_check< Mat<eT> > tmp2(A_tmp, eigvec);
00734 const Mat<eT>& A = tmp2.M;
00735
00736 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not square" );
00737
00738
00739
00740
00741
00742 char jobz = 'V';
00743 char uplo = 'U';
00744
00745 int n_rows = A.n_rows;
00746 int lwork = (std::max)(1, 3*n_rows-1);
00747
00748 eigval.set_size(n_rows);
00749 podarray<eT> work(lwork);
00750
00751 eigvec = A;
00752 int info;
00753
00754 arma_extra_debug_print("lapack::syev_()");
00755 lapack::syev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &n_rows, eigval.memptr(), work.memptr(), &lwork, &info);
00756 }
00757 #else
00758 {
00759 arma_stop("auxlib::eig_sym(): need LAPACK library");
00760 }
00761 #endif
00762
00763 }
00764
00765
00766
00767
00768 template<typename T>
00769 inline
00770 void
00771 auxlib::eig_sym(Col<T>& eigval, Mat< std::complex<T> >& eigvec, const Mat< std::complex<T> >& A_orig)
00772 {
00773 arma_extra_debug_sigprint();
00774
00775 typedef typename std::complex<T> eT;
00776
00777 #if defined(ARMA_USE_LAPACK)
00778 {
00779 const unwrap_check< Mat<eT> > tmp(A_orig, eigvec);
00780 const Mat<eT>& A = tmp.M;
00781
00782 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_sym(): given matrix is not hermitian" );
00783
00784 char jobz = 'V';
00785 char uplo = 'U';
00786
00787 int n_rows = A.n_rows;
00788 int lda = A.n_rows;
00789 int lwork = (std::max)(1, 2*n_rows - 1);
00790
00791 eigval.set_size(n_rows);
00792
00793 podarray<eT> work(lwork);
00794 podarray<T> rwork( (std::max)(1, 3*n_rows - 2) );
00795
00796 eigvec = A;
00797 int info;
00798
00799 arma_extra_debug_print("lapack::heev_()");
00800 lapack::heev_(&jobz, &uplo, &n_rows, eigvec.memptr(), &lda, eigval.memptr(), work.memptr(), &lwork, rwork.memptr(), &info);
00801 }
00802 #else
00803 {
00804 arma_stop("auxlib::eig_sym(): need LAPACK library");
00805 }
00806 #endif
00807
00808 }
00809
00810
00811
00812
00813
00814
00815 template<typename T>
00816 inline
00817 void
00818 auxlib::eig_gen
00819 (
00820 Col< std::complex<T> >& eigval,
00821 Mat<T>& l_eigvec,
00822 Mat<T>& r_eigvec,
00823 const Mat<T>& A,
00824 const char side
00825 )
00826 {
00827 arma_extra_debug_sigprint();
00828
00829
00830
00831 #if defined(ARMA_USE_LAPACK)
00832 {
00833 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" );
00834
00835 char jobvl;
00836 char jobvr;
00837
00838 switch(side)
00839 {
00840 case 'l':
00841 jobvl = 'V';
00842 jobvr = 'N';
00843 break;
00844
00845 case 'r':
00846 jobvl = 'N';
00847 jobvr = 'V';
00848 break;
00849
00850 case 'b':
00851 jobvl = 'V';
00852 jobvr = 'V';
00853 break;
00854
00855 case 'n':
00856 jobvl = 'N';
00857 jobvr = 'N';
00858 break;
00859
00860 default:
00861 arma_stop("auxlib::eig_gen(): parameter 'side' is invalid");
00862 }
00863
00864
00865 int n_rows = A.n_rows;
00866 int lda = A.n_rows;
00867 int lwork = (std::max)(1, 4*n_rows);
00868
00869 eigval.set_size(n_rows);
00870 l_eigvec.set_size(n_rows, n_rows);
00871 r_eigvec.set_size(n_rows, n_rows);
00872
00873 podarray<T> work(lwork);
00874 podarray<T> rwork( (std::max)(1, 3*n_rows) );
00875
00876 podarray<T> wr(n_rows);
00877 podarray<T> wi(n_rows);
00878
00879 Mat<T> A_copy = A;
00880 int info;
00881
00882 arma_extra_debug_print("lapack::cx_geev_()");
00883 lapack::geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, wr.memptr(), wi.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, &info);
00884
00885
00886 eigval.set_size(n_rows);
00887 for(u32 i=0; i<u32(n_rows); ++i)
00888 {
00889 eigval[i] = std::complex<T>(wr[i], wi[i]);
00890 }
00891
00892 }
00893 #else
00894 {
00895 arma_stop("auxlib::eig_gen(): need LAPACK library");
00896 }
00897 #endif
00898
00899 }
00900
00901
00902
00903
00904
00905
00906
00907
00908 template<typename T>
00909 inline
00910 void
00911 auxlib::eig_gen
00912 (
00913 Col< std::complex<T> >& eigval,
00914 Mat< std::complex<T> >& l_eigvec,
00915 Mat< std::complex<T> >& r_eigvec,
00916 const Mat< std::complex<T> >& A,
00917 const char side
00918 )
00919 {
00920 arma_extra_debug_sigprint();
00921
00922
00923
00924 typedef typename std::complex<T> eT;
00925
00926 #if defined(ARMA_USE_LAPACK)
00927 {
00928 arma_debug_check( (A.n_rows != A.n_cols), "auxlib::eig_gen(): given matrix is not square" );
00929
00930 char jobvl;
00931 char jobvr;
00932
00933 switch(side)
00934 {
00935 case 'l':
00936 jobvl = 'V';
00937 jobvr = 'N';
00938 break;
00939
00940 case 'r':
00941 jobvl = 'N';
00942 jobvr = 'V';
00943 break;
00944
00945 case 'b':
00946 jobvl = 'V';
00947 jobvr = 'V';
00948 break;
00949
00950 case 'n':
00951 jobvl = 'N';
00952 jobvr = 'N';
00953 break;
00954
00955 default:
00956 arma_stop("auxlib::eig_gen(): parameter 'side' is invalid");
00957 }
00958
00959
00960 int n_rows = A.n_rows;
00961 int lda = A.n_rows;
00962 int lwork = (std::max)(1, 4*n_rows);
00963
00964 eigval.set_size(n_rows);
00965 l_eigvec.set_size(n_rows, n_rows);
00966 r_eigvec.set_size(n_rows, n_rows);
00967
00968 podarray<eT> work(lwork);
00969 podarray<T> rwork( (std::max)(1, 3*n_rows) );
00970
00971 Mat<eT> A_copy = A;
00972 int info;
00973
00974 arma_extra_debug_print("lapack::cx_geev_()");
00975 lapack::cx_geev_(&jobvl, &jobvr, &n_rows, A_copy.memptr(), &lda, eigval.memptr(), l_eigvec.memptr(), &n_rows, r_eigvec.memptr(), &n_rows, work.memptr(), &lwork, rwork.memptr(), &info);
00976 }
00977 #else
00978 {
00979 arma_stop("auxlib::eig_gen(): need LAPACK library");
00980 }
00981 #endif
00982
00983 }
00984
00985
00986
00987 template<typename eT>
00988 inline
00989 bool
00990 auxlib::chol(Mat<eT>& out, const Mat<eT>& X)
00991 {
00992 arma_extra_debug_sigprint();
00993
00994
00995
00996
00997 #if defined(ARMA_USE_LAPACK)
00998 {
00999 char uplo = 'U';
01000 int n = X.n_rows;
01001 int lda = X.n_rows;
01002 int info;
01003
01004 out = X;
01005 lapack::potrf_(&uplo, &n, out.memptr(), &lda, &info);
01006
01007 for(u32 col=0; col<X.n_rows; ++col)
01008 {
01009 eT* colptr = out.colptr(col);
01010 for(u32 row=col+1; row<X.n_rows; ++row)
01011 {
01012 colptr[row] = eT(0);
01013 }
01014 }
01015
01016 return (info == 0);
01017 }
01018 #else
01019 {
01020 arma_stop("auxlib::chol(): need LAPACK library");
01021 return false;
01022 }
01023 #endif
01024 }
01025
01026
01027
01028 template<typename eT>
01029 inline
01030 bool
01031 auxlib::qr(Mat<eT>& Q, Mat<eT>& R, const Mat<eT>& X)
01032 {
01033 arma_extra_debug_sigprint();
01034
01035
01036
01037 #if defined(ARMA_USE_LAPACK)
01038 {
01039 int m = static_cast<int>(X.n_rows);
01040 int n = static_cast<int>(X.n_cols);
01041 int work_len = (std::max)(1,n);
01042 int work_len_tmp;
01043 int k = (std::min)(m,n);
01044 int info;
01045
01046 podarray<eT> tau(k);
01047 podarray<eT> work(work_len);
01048
01049 R = X;
01050
01051
01052 work_len_tmp = -1;
01053 lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01054
01055 if(info == 0)
01056 {
01057 work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01058 work.set_size(work_len);
01059 }
01060
01061 lapack::geqrf_(&m, &n, R.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01062
01063 Q.set_size(X.n_rows,X.n_rows);
01064
01065 eT* Q_mem = Q.memptr();
01066 const eT* R_mem = R.mem;
01067
01068 const u32 n_elem_copy = (std::min)(Q.n_elem, R.n_elem);
01069 for(u32 i=0; i < n_elem_copy; ++i)
01070 {
01071 Q_mem[i] = R_mem[i];
01072 }
01073
01074
01075
01076 for(u32 row=0; row < R.n_rows; ++row)
01077 {
01078 const u32 n_elem_tmp = (std::min)(row, R.n_cols);
01079 for(u32 col=0; col < n_elem_tmp; ++col)
01080 {
01081 R.at(row,col) = eT(0);
01082 }
01083 }
01084
01085
01086 if( (is_float<eT>::value == true) || (is_double<eT>::value == true) )
01087 {
01088
01089 work_len_tmp = -1;
01090 lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01091
01092 if(info == 0)
01093 {
01094 work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01095 work.set_size(work_len);
01096 }
01097
01098 lapack::orgqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01099 }
01100 else
01101 if( (is_supported_complex_float<eT>::value == true) || (is_supported_complex_double<eT>::value == true) )
01102 {
01103
01104 work_len_tmp = -1;
01105 lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len_tmp, &info);
01106
01107 if(info == 0)
01108 {
01109 work_len = static_cast<int>(auxlib::tmp_real(work[0]));
01110 work.set_size(work_len);
01111 }
01112
01113 lapack::ungqr_(&m, &m, &k, Q.memptr(), &m, tau.memptr(), work.memptr(), &work_len, &info);
01114 }
01115
01116 return (info == 0);
01117 }
01118 #else
01119 {
01120 arma_stop("auxlib::qr(): need LAPACK library");
01121 return false;
01122 }
01123 #endif
01124 }
01125
01126
01127
01128 template<typename eT>
01129 inline
01130 bool
01131 auxlib::svd(Col<eT>& S, const Mat<eT>& X)
01132 {
01133 arma_extra_debug_sigprint();
01134
01135
01136
01137 #if defined(ARMA_USE_LAPACK)
01138 {
01139 Mat<eT> A = X;
01140
01141 char jobu = 'N';
01142 char jobvt = 'N';
01143
01144 int m = A.n_rows;
01145 int n = A.n_cols;
01146 int lda = A.n_rows;
01147 int ldu = A.n_rows;
01148 int ldvt = A.n_cols;
01149 int lwork = 2;
01150 int info;
01151
01152 Mat<eT> U(1,1);
01153 Mat<eT> V(1,1);
01154
01155 S.set_size( (std::min)(m, n) );
01156
01157 podarray<eT> work(lwork);
01158
01159
01160
01161 int lwork_tmp = -1;
01162
01163 lapack::gesvd_<eT>
01164 (
01165 &jobu, &jobvt,
01166 &m,&n,
01167 A.memptr(), &lda,
01168 S.memptr(),
01169 U.memptr(), &ldu,
01170 V.memptr(), &ldvt,
01171 work.memptr(), &lwork_tmp,
01172 &info
01173 );
01174
01175 if(info == 0)
01176 {
01177 lwork = static_cast<int>(work[0]);
01178 work.set_size(lwork);
01179
01180 lapack::gesvd_<eT>
01181 (
01182 &jobu, &jobvt,
01183 &m, &n,
01184 A.memptr(), &lda,
01185 S.memptr(),
01186 U.memptr(), &ldu,
01187 V.memptr(), &ldvt,
01188 work.memptr(), &lwork,
01189 &info
01190 );
01191
01192 return (info == 0);
01193 }
01194 else
01195 {
01196 return false;
01197 }
01198 }
01199 #else
01200 {
01201 arma_stop("auxlib::svd(): need LAPACK library");
01202 return false;
01203 }
01204 #endif
01205 }
01206
01207
01208
01209 template<typename T>
01210 inline
01211 bool
01212 auxlib::svd(Col<T>& S, const Mat< std::complex<T> >& X)
01213 {
01214 arma_extra_debug_sigprint();
01215
01216 typedef std::complex<T> eT;
01217
01218
01219
01220 #if defined(ARMA_USE_LAPACK)
01221 {
01222 Mat<eT> A = X;
01223
01224 char jobu = 'N';
01225 char jobvt = 'N';
01226
01227 int m = A.n_rows;
01228 int n = A.n_cols;
01229 int lda = A.n_rows;
01230 int ldu = A.n_rows;
01231 int ldvt = A.n_cols;
01232 int lwork = 2 * (std::min)(m,n) + (std::max)(m,n);
01233 int info;
01234
01235 Mat<eT> U(1,1);
01236 Mat<eT> V(1,1);
01237
01238 S.set_size( (std::min)(m,n) );
01239
01240 podarray<eT> work(lwork);
01241 podarray<T> rwork( 5*(std::min)(m,n) );
01242
01243
01244 int lwork_tmp = -1;
01245
01246 lapack::cx_gesvd_<T>
01247 (
01248 &jobu, &jobvt,
01249 &m, &n,
01250 A.memptr(), &lda,
01251 S.memptr(),
01252 U.memptr(), &ldu,
01253 V.memptr(), &ldvt,
01254 work.memptr(), &lwork_tmp,
01255 rwork.memptr(),
01256 &info
01257 );
01258
01259 if(info == 0)
01260 {
01261 int proposed_lwork = static_cast<int>(real(work[0]));
01262 if(proposed_lwork > lwork)
01263 {
01264 lwork = proposed_lwork;
01265 work.set_size(lwork);
01266 }
01267
01268 lapack::cx_gesvd_<T>
01269 (
01270 &jobu, &jobvt,
01271 &m, &n,
01272 A.memptr(), &lda,
01273 S.memptr(),
01274 U.memptr(), &ldu,
01275 V.memptr(), &ldvt,
01276 work.memptr(), &lwork,
01277 rwork.memptr(),
01278 &info
01279 );
01280
01281 return (info == 0);
01282 }
01283 else
01284 {
01285 return false;
01286 }
01287 }
01288 #else
01289 {
01290 arma_stop("auxlib::svd(): need LAPACK library");
01291 return false;
01292 }
01293 #endif
01294 }
01295
01296
01297
01298 template<typename eT>
01299 inline
01300 bool
01301 auxlib::svd(Mat<eT>& U, Col<eT>& S, Mat<eT>& V, const Mat<eT>& X)
01302 {
01303 arma_extra_debug_sigprint();
01304
01305
01306
01307 #if defined(ARMA_USE_LAPACK)
01308 {
01309 Mat<eT> A = X;
01310
01311 char jobu = 'A';
01312 char jobvt = 'A';
01313
01314 int m = A.n_rows;
01315 int n = A.n_cols;
01316 int lda = A.n_rows;
01317 int ldu = A.n_rows;
01318 int ldvt = A.n_cols;
01319 int lwork = 2;
01320 int info;
01321
01322 U.set_size(m,m);
01323 V.set_size(n,n);
01324
01325 S.set_size( (std::min)(m,n) );
01326 podarray<eT> work(lwork);
01327
01328
01329 int lwork_tmp = -1;
01330
01331 lapack::gesvd_<eT>
01332 (
01333 &jobu, &jobvt,
01334 &m, &n,
01335 A.memptr(), &lda,
01336 S.memptr(),
01337 U.memptr(), &ldu,
01338 V.memptr(), &ldvt,
01339 work.memptr(), &lwork_tmp,
01340 &info
01341 );
01342
01343 if(info == 0)
01344 {
01345 lwork = static_cast<int>(work[0]);
01346 work.set_size(lwork);
01347
01348 lapack::gesvd_<eT>
01349 (
01350 &jobu, &jobvt,
01351 &m, &n,
01352 A.memptr(), &lda,
01353 S.memptr(),
01354 U.memptr(), &ldu,
01355 V.memptr(), &ldvt,
01356 work.memptr(), &lwork,
01357 &info
01358 );
01359
01360 op_trans::apply(V,V);
01361
01362 return (info == 0);
01363 }
01364 else
01365 {
01366 return false;
01367 }
01368 }
01369 #else
01370 {
01371 arma_stop("auxlib::svd(): need LAPACK library");
01372 return false;
01373 }
01374 #endif
01375 }
01376
01377
01378
01379 template<typename T>
01380 inline
01381 bool
01382 auxlib::svd(Mat< std::complex<T> >& U, Col<T> &S, Mat< std::complex<T> >& V, const Mat< std::complex<T> >& X)
01383 {
01384 arma_extra_debug_sigprint();
01385
01386 typedef std::complex<T> eT;
01387
01388
01389
01390 #if defined(ARMA_USE_LAPACK)
01391 {
01392 Mat<eT> A = X;
01393
01394 char jobu = 'A';
01395 char jobvt = 'A';
01396
01397 int m = A.n_rows;
01398 int n = A.n_cols;
01399 int lda = A.n_rows;
01400 int ldu = A.n_rows;
01401 int ldvt = A.n_cols;
01402 int lwork = 2;
01403 int info;
01404
01405 U.set_size(m,m);
01406 V.set_size(n,n);
01407
01408 S.set_size( (std::min)(m,n) );
01409
01410 podarray<eT> work(lwork);
01411 podarray<T> rwork( 5*(std::min)(m,n) );
01412
01413
01414 int lwork_tmp = -1;
01415 lapack::cx_gesvd_<T>
01416 (
01417 &jobu, &jobvt,
01418 &m, &n,
01419 A.memptr(), &lda,
01420 S.memptr(),
01421 U.memptr(), &ldu,
01422 V.memptr(), &ldvt,
01423 work.memptr(), &lwork_tmp,
01424 rwork.memptr(),
01425 &info
01426 );
01427
01428 if(info == 0)
01429 {
01430 lwork = static_cast<int>(real(work[0]));
01431 work.set_size(lwork);
01432
01433 lapack::cx_gesvd_<T>
01434 (
01435 &jobu, &jobvt,
01436 &m, &n,
01437 A.memptr(), &lda,
01438 S.memptr(),
01439 U.memptr(), &ldu,
01440 V.memptr(), &ldvt,
01441 work.memptr(), &lwork,
01442 rwork.memptr(),
01443 &info
01444 );
01445
01446 op_htrans::apply(V,V);
01447
01448 return (info == 0);
01449 }
01450 else
01451 {
01452 return false;
01453 }
01454 }
01455 #else
01456 {
01457 arma_stop("auxlib::svd(): need LAPACK library");
01458 return false;
01459 }
01460 #endif
01461
01462 }
01463
01464
01465
01466
01467
01468 template<typename eT>
01469 inline
01470 bool
01471 auxlib::solve(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
01472 {
01473 arma_extra_debug_sigprint();
01474
01475 #if defined(ARMA_USE_LAPACK)
01476 {
01477 int n = A.n_rows;
01478 int lda = A.n_rows;
01479 int ldb = A.n_rows;
01480 int nrhs = B.n_cols;
01481 int info;
01482
01483 podarray<int> ipiv(n);
01484
01485 out = B;
01486 Mat<eT> A_copy = A;
01487
01488 lapack::gesv_<eT>(&n, &nrhs, A_copy.memptr(), &lda, ipiv.memptr(), out.memptr(), &ldb, &info);
01489
01490 return (info == 0);
01491 }
01492 #else
01493 {
01494 arma_stop("auxlib::solve(): need LAPACK library");
01495 return false;
01496 }
01497 #endif
01498 }
01499
01500
01501
01502
01503
01504
01505 template<typename eT>
01506 inline
01507 bool
01508 auxlib::solve_od(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
01509 {
01510 arma_extra_debug_sigprint();
01511
01512 #if defined(ARMA_USE_LAPACK)
01513 {
01514 char trans = 'N';
01515
01516 int m = A.n_rows;
01517 int n = A.n_cols;
01518 int lda = A.n_rows;
01519 int ldb = A.n_rows;
01520 int nrhs = B.n_cols;
01521 int lwork = n + (std::max)(n, nrhs);
01522 int info;
01523
01524 Mat<eT> A_copy = A;
01525 Mat<eT> tmp = B;
01526
01527
01528 podarray<eT> work(lwork);
01529
01530 arma_extra_debug_print("lapack::gels_()");
01531
01532
01533
01534
01535
01536 lapack::gels_<eT>
01537 (
01538 &trans, &m, &n, &nrhs,
01539 A_copy.memptr(), &lda,
01540 tmp.memptr(), &ldb,
01541 work.memptr(), &lwork,
01542 &info
01543 );
01544
01545 arma_extra_debug_print("lapack::gels_() -- finished");
01546
01547 out.set_size(A.n_cols, B.n_cols);
01548
01549 for(u32 col=0; col<B.n_cols; ++col)
01550 {
01551 syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols );
01552 }
01553
01554 return (info == 0);
01555 }
01556 #else
01557 {
01558 arma_stop("auxlib::solve_od(): need LAPACK library");
01559 return false;
01560 }
01561 #endif
01562 }
01563
01564
01565
01566
01567
01568
01569 template<typename eT>
01570 inline
01571 bool
01572 auxlib::solve_ud(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B)
01573 {
01574 arma_extra_debug_sigprint();
01575
01576 #if defined(ARMA_USE_LAPACK)
01577 {
01578 char trans = 'N';
01579
01580 int m = A.n_rows;
01581 int n = A.n_cols;
01582 int lda = A.n_rows;
01583 int ldb = A.n_cols;
01584 int nrhs = B.n_cols;
01585 int lwork = m + (std::max)(m,nrhs);
01586 int info;
01587
01588
01589 Mat<eT> A_copy = A;
01590
01591 Mat<eT> tmp;
01592 tmp.zeros(A.n_cols, B.n_cols);
01593
01594 for(u32 col=0; col<B.n_cols; ++col)
01595 {
01596 eT* tmp_colmem = tmp.colptr(col);
01597
01598 syslib::copy_elem( tmp_colmem, B.colptr(col), B.n_rows );
01599
01600 for(u32 row=B.n_rows; row<A.n_cols; ++row)
01601 {
01602 tmp_colmem[row] = eT(0);
01603 }
01604 }
01605
01606 podarray<eT> work(lwork);
01607
01608 arma_extra_debug_print("lapack::gels_()");
01609
01610
01611
01612
01613
01614 lapack::gels_<eT>
01615 (
01616 &trans, &m, &n, &nrhs,
01617 A_copy.memptr(), &lda,
01618 tmp.memptr(), &ldb,
01619 work.memptr(), &lwork,
01620 &info
01621 );
01622
01623 arma_extra_debug_print("lapack::gels_() -- finished");
01624
01625 out.set_size(A.n_cols, B.n_cols);
01626
01627 for(u32 col=0; col<B.n_cols; ++col)
01628 {
01629 syslib::copy_elem( out.colptr(col), tmp.colptr(col), A.n_cols );
01630 }
01631
01632 return (info == 0);
01633 }
01634 #else
01635 {
01636 arma_stop("auxlib::solve_ud(): need LAPACK library");
01637 return false;
01638 }
01639 #endif
01640 }
01641
01642
01643
01644