00001 00033 #include <itpp/base/matfunc.h> 00034 #include <itpp/base/schur.h> 00035 #include <itpp/base/converters.h> 00036 #include <limits> 00037 00038 00039 namespace itpp { 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 // ---------------------- Instantiations ------------------------------------ 00082 00083 template int length(const vec &v); 00084 template int length(const cvec &v); 00085 template int length(const svec &v); 00086 template int length(const ivec &v); 00087 template int length(const bvec &v); 00088 00089 00090 template double sum(const vec &v); 00091 template std::complex<double> sum(const cvec &v); 00092 template short sum(const svec &v); 00093 template int sum(const ivec &v); 00094 template bin sum(const bvec &v); 00095 00096 template double sum_sqr(const vec &v); 00097 template std::complex<double> sum_sqr(const cvec &v); 00098 template short sum_sqr(const svec &v); 00099 template int sum_sqr(const ivec &v); 00100 template bin sum_sqr(const bvec &v); 00101 00102 template vec cumsum(const vec &v); 00103 template cvec cumsum(const cvec &v); 00104 template svec cumsum(const svec &v); 00105 template ivec cumsum(const ivec &v); 00106 template bvec cumsum(const bvec &v); 00107 00108 template double prod(const vec &v); 00109 template std::complex<double> prod(const cvec &v); 00110 template short prod(const svec &v); 00111 template int prod(const ivec &v); 00112 template bin prod(const bvec &v); 00113 00114 template vec cross(const vec &v1, const vec &v2); 00115 template ivec cross(const ivec &v1, const ivec &v2); 00116 template svec cross(const svec &v1, const svec &v2); 00117 00118 template vec reverse(const vec &in); 00119 template cvec reverse(const cvec &in); 00120 template svec reverse(const svec &in); 00121 template ivec reverse(const ivec &in); 00122 template bvec reverse(const bvec &in); 00123 00124 template vec repeat(const vec &v, int norepeats); 00125 template cvec repeat(const cvec &v, int norepeats); 00126 template svec repeat(const svec &v, int norepeats); 00127 template ivec repeat(const ivec &v, int norepeats); 00128 template bvec repeat(const bvec &v, int norepeats); 00129 00130 template vec apply_function(float (*f)(float), const vec &data); 00131 template vec apply_function(double (*f)(double), const vec &data); 00132 template cvec apply_function(std::complex<double> (*f)(std::complex<double>), const cvec &data); 00133 template svec apply_function(short (*f)(short), const svec &data); 00134 template ivec apply_function(int (*f)(int), const ivec &data); 00135 template bvec apply_function(bin (*f)(bin), const bvec &data); 00136 00137 00138 template ivec zero_pad(const ivec &v, int n); 00139 template vec zero_pad(const vec &v, int n); 00140 template cvec zero_pad(const cvec &v, int n); 00141 template bvec zero_pad(const bvec &v, int n); 00142 00143 template ivec zero_pad(const ivec &v); 00144 template vec zero_pad(const vec &v); 00145 template cvec zero_pad(const cvec &v); 00146 template bvec zero_pad(const bvec &v); 00147 00148 template mat zero_pad(const mat &, int, int); 00149 template cmat zero_pad(const cmat &, int, int); 00150 template imat zero_pad(const imat &, int, int); 00151 template bmat zero_pad(const bmat &, int, int); 00152 00153 template vec sum(const mat &m, int dim); 00154 template cvec sum(const cmat &m, int dim); 00155 template svec sum(const smat &m, int dim); 00156 template ivec sum(const imat &m, int dim); 00157 template bvec sum(const bmat &m, int dim); 00158 00159 template vec sum_sqr(const mat & m, int dim); 00160 template cvec sum_sqr(const cmat &m, int dim); 00161 template svec sum_sqr(const smat &m, int dim); 00162 template ivec sum_sqr(const imat &m, int dim); 00163 template bvec sum_sqr(const bmat &m, int dim); 00164 00165 template mat cumsum(const mat &m, int dim); 00166 template cmat cumsum(const cmat &m, int dim); 00167 template smat cumsum(const smat &m, int dim); 00168 template imat cumsum(const imat &m, int dim); 00169 template bmat cumsum(const bmat &m, int dim); 00170 00171 template vec prod(const mat &m, int dim); 00172 // Template instantiation of product 00173 template cvec prod(const cmat &v, int dim); 00174 template svec prod(const smat &m, int dim); 00175 template ivec prod(const imat &m, int dim); 00176 00177 template vec diag(const mat &in); 00178 template cvec diag(const cmat &in); 00179 00180 template void diag(const vec &in, mat &m); 00181 template void diag(const cvec &in, cmat &m); 00182 00183 template mat diag(const vec &v, const int K); 00184 template cmat diag(const cvec &v, const int K); 00185 00186 template mat bidiag(const vec &, const vec &); 00187 template cmat bidiag(const cvec &, const cvec &); 00188 00189 template void bidiag(const vec &, const vec &, mat &); 00190 template void bidiag(const cvec &, const cvec &, cmat &); 00191 00192 template void bidiag(const mat &, vec &, vec &); 00193 template void bidiag(const cmat &, cvec &, cvec &); 00194 00195 template mat tridiag(const vec &main, const vec &, const vec &); 00196 template cmat tridiag(const cvec &main, const cvec &, const cvec &); 00197 00198 template void tridiag(const vec &main, const vec &, const vec &, mat &); 00199 template void tridiag(const cvec &main, const cvec &, const cvec &, cmat &); 00200 00201 template void tridiag(const mat &m, vec &, vec &, vec &); 00202 template void tridiag(const cmat &m, cvec &, cvec &, cvec &); 00203 00204 template double trace(const mat &in); 00205 template std::complex<double> trace(const cmat &in); 00206 template short trace(const smat &in); 00207 template int trace(const imat &in); 00208 template bin trace(const bmat &in); 00209 00210 template void transpose(const mat &m, mat &out); 00211 template void transpose(const cmat &m, cmat &out); 00212 template void transpose(const smat &m, smat &out); 00213 template void transpose(const imat &m, imat &out); 00214 template void transpose(const bmat &m, bmat &out); 00215 00216 template mat transpose(const mat &m); 00217 template cmat transpose(const cmat &m); 00218 template smat transpose(const smat &m); 00219 template imat transpose(const imat &m); 00220 template bmat transpose(const bmat &m); 00221 00222 00223 template void hermitian_transpose(const mat &m, mat &out); 00224 template void hermitian_transpose(const cmat &m, cmat &out); 00225 template void hermitian_transpose(const smat &m, smat &out); 00226 template void hermitian_transpose(const imat &m, imat &out); 00227 template void hermitian_transpose(const bmat &m, bmat &out); 00228 00229 template mat hermitian_transpose(const mat &m); 00230 template cmat hermitian_transpose(const cmat &m); 00231 template smat hermitian_transpose(const smat &m); 00232 template imat hermitian_transpose(const imat &m); 00233 template bmat hermitian_transpose(const bmat &m); 00234 00235 template mat repeat(const mat &m, int norepeats); 00236 template cmat repeat(const cmat &m, int norepeats); 00237 template smat repeat(const smat &m, int norepeats); 00238 template imat repeat(const imat &m, int norepeats); 00239 template bmat repeat(const bmat &m, int norepeats); 00240 00241 template mat apply_function(float (*f)(float), const mat &data); 00242 template mat apply_function(double (*f)(double), const mat &data); 00243 template cmat apply_function(std::complex<double> (*f)(std::complex<double>), const cmat &data); 00244 template smat apply_function(short (*f)(short), const smat &data); 00245 template imat apply_function(int (*f)(int), const imat &data); 00246 template bmat apply_function(bin (*f)(bin), const bmat &data); 00247 00248 template vec rvectorize(const mat &m); 00249 template cvec rvectorize(const cmat &m); 00250 template ivec rvectorize(const imat &m); 00251 template bvec rvectorize(const bmat &m); 00252 00253 template vec cvectorize(const mat &m); 00254 template cvec cvectorize(const cmat &m); 00255 template ivec cvectorize(const imat &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 bmat reshape(const bmat &m, int rows, int cols); 00262 00263 template mat reshape(const vec &m, int rows, int cols); 00264 template cmat reshape(const cvec &m, int rows, int cols); 00265 template imat reshape(const ivec &m, int rows, int cols); 00266 template bmat reshape(const bvec &m, int rows, int cols); 00267 00268 template vec upsample(const vec &v, int usf); 00269 template cvec upsample(const cvec &v, int usf); 00270 template svec upsample(const svec &v, int usf); 00271 template ivec upsample(const ivec &v, int usf); 00272 template bvec upsample(const bvec &v, int usf); 00273 00274 template mat upsample(const mat &v, int usf); 00275 template cmat upsample(const cmat &v, int usf); 00276 template smat upsample(const smat &v, int usf); 00277 template imat upsample(const imat &v, int usf); 00278 template bmat upsample(const bmat &v, int usf); 00279 00280 template void upsample(const vec &v, int usf, vec & u); 00281 template void upsample(const cvec &v, int usf, cvec & u); 00282 template void upsample(const svec &v, int usf, svec & u); 00283 template void upsample(const ivec &v, int usf, ivec & u); 00284 template void upsample(const bvec &v, int usf, bvec & u); 00285 00286 template void upsample(const mat &v, int usf, mat & u); 00287 template void upsample(const cmat &v, int usf, cmat & u); 00288 template void upsample(const smat &v, int usf, smat & u); 00289 template void upsample(const imat &v, int usf, imat & u); 00290 template void upsample(const bmat &v, int usf, bmat & u); 00291 00292 template vec lininterp(const vec &v, int usf); 00293 template cvec lininterp(const cvec &v, int usf); 00294 00295 template mat lininterp(const mat &v, int usf); 00296 template cmat lininterp(const cmat &v, int usf); 00297 00298 template void lininterp(const vec &v, int usf, vec & u); 00299 template void lininterp(const cvec &v, int usf, cvec & u); 00300 00301 template void lininterp(const mat &v, int usf, mat & u); 00302 template void lininterp(const cmat &v, int usf, cmat & u); 00303 00304 template mat lininterp(const mat &m, const double f_base, const double f_ups, const int nrof_samples, const double t_start); 00305 template cmat lininterp(const cmat &m, const double f_base, const double f_ups, const int nrof_samples, const double t_start); 00306 00307 template vec lininterp(const vec &v, const double f_base, const double f_ups, const int nrof_samples, const double t_start); 00308 template cvec lininterp(const cvec &v, const double f_base, const double f_ups, const int nrof_samples, const double t_start); 00309 00310 } // namespace itpp
Generated on Wed Apr 18 11:23:30 2007 for IT++ by Doxygen 1.5.2