00001 00031 #include <itpp/base/matfunc.h> 00032 #include <itpp/base/algebra/schur.h> 00033 #include <itpp/base/converters.h> 00034 #include <limits> 00035 00037 00038 namespace itpp 00039 { 00040 00041 // Square root of a square matrix (based on Octave sqrtm implementation) 00042 cmat sqrtm(const mat& A) 00043 { 00044 return sqrtm(to_cmat(A)); 00045 } 00046 00047 // Square root of the complex square matrix A 00048 cmat sqrtm(const cmat& A) 00049 { 00050 cmat U, T; 00051 schur(A, U, T); 00052 00053 int n = U.rows(); 00054 cmat R(n, n); 00055 00056 R.zeros(); 00057 for (int j = 0; j < n; j++) 00058 R(j, j) = std::sqrt(T(j, j)); 00059 00060 const double fudge = std::sqrt(std::numeric_limits<double>::min()); 00061 00062 for (int p = 0; p < n - 1; p++) { 00063 for (int i = 0; i < n - (p + 1); i++) { 00064 const int j = i + p + 1; 00065 std::complex<double> s = T(i, j); 00066 00067 for (int k = i + 1; k < j; k++) 00068 s -= R(i, k) * R(k, j); 00069 00070 const std::complex<double> d = R(i, i) + R(j, j) + fudge; 00071 const std::complex<double> conj_d = conj(d); 00072 00073 R(i, j) = (s * conj_d) / (d * conj_d); 00074 } 00075 } 00076 00077 return U * R * U.H(); 00078 } 00079 00080 00081 bool all(const bvec &testvec) 00082 { 00083 for (int i = 0; i < testvec.length(); i++) 00084 if (!testvec(i)) return false; 00085 return true; 00086 } 00087 00088 bool any(const bvec &testvec) 00089 { 00090 for (int i = 0; i < testvec.length(); i++) 00091 if (testvec(i)) return true; 00092 return false; 00093 } 00094 00095 // ---------------------------------------------------------------------- 00096 // Instantiations 00097 // ---------------------------------------------------------------------- 00098 00099 template int length(const vec &v); 00100 template int length(const cvec &v); 00101 template int length(const svec &v); 00102 template int length(const ivec &v); 00103 template int length(const bvec &v); 00104 00105 template double sum(const vec &v); 00106 template std::complex<double> sum(const cvec &v); 00107 template short sum(const svec &v); 00108 template int sum(const ivec &v); 00109 template bin sum(const bvec &v); 00110 00111 template double sum_sqr(const vec &v); 00112 template std::complex<double> sum_sqr(const cvec &v); 00113 template short sum_sqr(const svec &v); 00114 template int sum_sqr(const ivec &v); 00115 template bin sum_sqr(const bvec &v); 00116 00117 template vec cumsum(const vec &v); 00118 template cvec cumsum(const cvec &v); 00119 template svec cumsum(const svec &v); 00120 template ivec cumsum(const ivec &v); 00121 template bvec cumsum(const bvec &v); 00122 00123 template double prod(const vec &v); 00124 template std::complex<double> prod(const cvec &v); 00125 template short prod(const svec &v); 00126 template int prod(const ivec &v); 00127 template bin prod(const bvec &v); 00128 00129 template vec cross(const vec &v1, const vec &v2); 00130 template cvec cross(const cvec &v1, const cvec &v2); 00131 template ivec cross(const ivec &v1, const ivec &v2); 00132 template svec cross(const svec &v1, const svec &v2); 00133 template bvec cross(const bvec &v1, const bvec &v2); 00134 00135 template vec reverse(const vec &in); 00136 template cvec reverse(const cvec &in); 00137 template svec reverse(const svec &in); 00138 template ivec reverse(const ivec &in); 00139 template bvec reverse(const bvec &in); 00140 00141 template vec zero_pad(const vec &v, int n); 00142 template cvec zero_pad(const cvec &v, int n); 00143 template ivec zero_pad(const ivec &v, int n); 00144 template svec zero_pad(const svec &v, int n); 00145 template bvec zero_pad(const bvec &v, int n); 00146 00147 template vec zero_pad(const vec &v); 00148 template cvec zero_pad(const cvec &v); 00149 template ivec zero_pad(const ivec &v); 00150 template svec zero_pad(const svec &v); 00151 template bvec zero_pad(const bvec &v); 00152 00153 template mat zero_pad(const mat &, int, int); 00154 template cmat zero_pad(const cmat &, int, int); 00155 template imat zero_pad(const imat &, int, int); 00156 template smat zero_pad(const smat &, int, int); 00157 template bmat zero_pad(const bmat &, int, int); 00158 00159 template vec sum(const mat &m, int dim); 00160 template cvec sum(const cmat &m, int dim); 00161 template svec sum(const smat &m, int dim); 00162 template ivec sum(const imat &m, int dim); 00163 template bvec sum(const bmat &m, int dim); 00164 00165 template double sumsum(const mat &X); 00166 template std::complex<double> sumsum(const cmat &X); 00167 template short sumsum(const smat &X); 00168 template int sumsum(const imat &X); 00169 template bin sumsum(const bmat &X); 00170 00171 template vec sum_sqr(const mat & m, int dim); 00172 template cvec sum_sqr(const cmat &m, int dim); 00173 template svec sum_sqr(const smat &m, int dim); 00174 template ivec sum_sqr(const imat &m, int dim); 00175 template bvec sum_sqr(const bmat &m, int dim); 00176 00177 template mat cumsum(const mat &m, int dim); 00178 template cmat cumsum(const cmat &m, int dim); 00179 template smat cumsum(const smat &m, int dim); 00180 template imat cumsum(const imat &m, int dim); 00181 template bmat cumsum(const bmat &m, int dim); 00182 00183 template vec prod(const mat &m, int dim); 00184 template cvec prod(const cmat &v, int dim); 00185 template svec prod(const smat &m, int dim); 00186 template ivec prod(const imat &m, int dim); 00187 template bvec prod(const bmat &m, int dim); 00188 00189 template vec diag(const mat &in); 00190 template cvec diag(const cmat &in); 00191 00192 template void diag(const vec &in, mat &m); 00193 template void diag(const cvec &in, cmat &m); 00194 00195 template mat diag(const vec &v, const int K); 00196 template cmat diag(const cvec &v, const int K); 00197 00198 template mat bidiag(const vec &, const vec &); 00199 template cmat bidiag(const cvec &, const cvec &); 00200 00201 template void bidiag(const vec &, const vec &, mat &); 00202 template void bidiag(const cvec &, const cvec &, cmat &); 00203 00204 template void bidiag(const mat &, vec &, vec &); 00205 template void bidiag(const cmat &, cvec &, cvec &); 00206 00207 template mat tridiag(const vec &main, const vec &, const vec &); 00208 template cmat tridiag(const cvec &main, const cvec &, const cvec &); 00209 00210 template void tridiag(const vec &main, const vec &, const vec &, mat &); 00211 template void tridiag(const cvec &main, const cvec &, const cvec &, cmat &); 00212 00213 template void tridiag(const mat &m, vec &, vec &, vec &); 00214 template void tridiag(const cmat &m, cvec &, cvec &, cvec &); 00215 00216 template double trace(const mat &in); 00217 template std::complex<double> trace(const cmat &in); 00218 template short trace(const smat &in); 00219 template int trace(const imat &in); 00220 template bin trace(const bmat &in); 00221 00222 template void transpose(const mat &m, mat &out); 00223 template void transpose(const cmat &m, cmat &out); 00224 template void transpose(const smat &m, smat &out); 00225 template void transpose(const imat &m, imat &out); 00226 template void transpose(const bmat &m, bmat &out); 00227 00228 template mat transpose(const mat &m); 00229 template cmat transpose(const cmat &m); 00230 template smat transpose(const smat &m); 00231 template imat transpose(const imat &m); 00232 template bmat transpose(const bmat &m); 00233 00234 template void hermitian_transpose(const mat &m, mat &out); 00235 template void hermitian_transpose(const cmat &m, cmat &out); 00236 template void hermitian_transpose(const smat &m, smat &out); 00237 template void hermitian_transpose(const imat &m, imat &out); 00238 template void hermitian_transpose(const bmat &m, bmat &out); 00239 00240 template mat hermitian_transpose(const mat &m); 00241 template cmat hermitian_transpose(const cmat &m); 00242 template smat hermitian_transpose(const smat &m); 00243 template imat hermitian_transpose(const imat &m); 00244 template bmat hermitian_transpose(const bmat &m); 00245 00246 template vec rvectorize(const mat &m); 00247 template cvec rvectorize(const cmat &m); 00248 template ivec rvectorize(const imat &m); 00249 template svec rvectorize(const smat &m); 00250 template bvec rvectorize(const bmat &m); 00251 00252 template vec cvectorize(const mat &m); 00253 template cvec cvectorize(const cmat &m); 00254 template ivec cvectorize(const imat &m); 00255 template svec cvectorize(const smat &m); 00256 template bvec cvectorize(const bmat &m); 00257 00258 template mat reshape(const mat &m, int rows, int cols); 00259 template cmat reshape(const cmat &m, int rows, int cols); 00260 template imat reshape(const imat &m, int rows, int cols); 00261 template smat reshape(const smat &m, int rows, int cols); 00262 template bmat reshape(const bmat &m, int rows, int cols); 00263 00264 template mat reshape(const vec &m, int rows, int cols); 00265 template cmat reshape(const cvec &m, int rows, int cols); 00266 template imat reshape(const ivec &m, int rows, int cols); 00267 template smat reshape(const svec &m, int rows, int cols); 00268 template bmat reshape(const bvec &m, int rows, int cols); 00269 00270 template mat kron(const mat &X, const mat &Y); 00271 template cmat kron(const cmat &X, const cmat &Y); 00272 template imat kron(const imat &X, const imat &Y); 00273 template smat kron(const smat &X, const smat &Y); 00274 template bmat kron(const bmat &X, const bmat &Y); 00275 00276 template vec repmat(const vec &v, int n); 00277 template cvec repmat(const cvec &v, int n); 00278 template ivec repmat(const ivec &v, int n); 00279 template svec repmat(const svec &v, int n); 00280 template bvec repmat(const bvec &v, int n); 00281 00282 template mat repmat(const vec &v, int m, int n, bool transpose); 00283 template cmat repmat(const cvec &v, int m, int n, bool transpose); 00284 template imat repmat(const ivec &v, int m, int n, bool transpose); 00285 template smat repmat(const svec &v, int m, int n, bool transpose); 00286 template bmat repmat(const bvec &v, int m, int n, bool transpose); 00287 00288 template mat repmat(const mat &data, int m, int n); 00289 template cmat repmat(const cmat &data, int m, int n); 00290 template imat repmat(const imat &data, int m, int n); 00291 template smat repmat(const smat &data, int m, int n); 00292 template bmat repmat(const bmat &data, int m, int n); 00293 00294 } // namespace itpp 00295
Generated on Wed Mar 2 2011 22:05:07 for IT++ by Doxygen 1.7.3