Lapack++
|
00001 /*-*-c++-*-**************************************************************** 00002 * latmpl.h C++ Templates for lapackpp classes 00003 ------------------- 00004 begin : 2005-12-29 00005 copyright : (C) 2005 by Christian Stimming 00006 email : stimming@tuhh.de 00007 ***************************************************************************/ 00008 00009 // This library is free software; you can redistribute it and/or 00010 // modify it under the terms of the GNU Lesser General Public License as 00011 // published by the Free Software Foundation; either version 2, or (at 00012 // your option) any later version. 00013 00014 // This library is distributed in the hope that it will be useful, 00015 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00017 // GNU Lesser General Public License for more details. 00018 00019 // You should have received a copy of the GNU Lesser General Public License along 00020 // with this library; see the file COPYING. If not, write to the Free 00021 // Software Foundation, 59 Temple Place - Suite 330, Boston, MA 02111-1307, 00022 // USA. 00023 00024 #ifndef _LATMPL_H 00025 #define _LATMPL_H 00026 00034 #include <cstdlib> 00035 #include <cmath> 00036 #include <lafnames.h> 00037 #include LA_EXCEPTION_H 00038 00039 // Watch out: Only some compilers correctly implement the C++ 00040 // standard specifications, where type names from inside template 00041 // classes need to be preceded by the "typename" keyword. Other 00042 // compilers do not accept the keyword "typename". For those we 00043 // define it as empty here. 00044 #ifndef __GNUC__ 00045 // Non-gcc compiler 00046 # ifdef _MSC_VER 00047 // This is Microsoft Visual C++. 00048 // Microsoft Visual C++ 7.1 _MSC_VER = 1310 00049 // Microsoft Visual C++ 7.0 _MSC_VER = 1300 00050 // Microsoft Visual C++ 6.0 _MSC_VER = 1200 00051 // Microsoft Visual C++ 5.0 _MSC_VER = 1100 00052 # if _MSC_VER <= 1300 00053 // This is Microsoft Visual C++ 7.0 or older, so define the 00054 // keyword as empty 00055 # define typename 00056 # else 00057 // Microsoft Visual C++ 7.1 or newer. typename should be ok. 00058 # endif 00059 # else 00060 // Compiler unknown 00061 # endif 00062 #endif 00063 00064 namespace la 00065 { 00066 00070 template <class matT> 00071 bool is_zero(const matT& mat) 00072 { 00073 int i, j, M = mat.rows(), N = mat.cols(); 00074 00075 // If your compiler gives an error in this line, please see the 00076 // note above on "typename". 00077 typename matT::value_type zero(0); 00078 00079 for (j = 0; j < N; j++) 00080 for (i = 0; i < M; i++) 00081 if (mat(i, j) != zero) 00082 return false; 00083 return true; 00084 } 00085 00088 template <class matT> 00089 bool equal(const matT& mat1, const matT& mat2) 00090 { 00091 int i, j, M = mat1.rows(), N = mat1.cols(); 00092 if (mat1.rows() != mat2.rows() || mat1.cols() != mat2.cols()) 00093 throw LaException("equal<matT>(const matT&, const matT)", 00094 "Both matrices have unequal sizes"); 00095 for (j = 0; j < N; j++) 00096 for (i = 0; i < M; i++) 00097 if (mat1(i, j) != mat2(i, j)) 00098 return false; 00099 return true; 00100 } 00102 00103 00109 template <class matT> 00110 void zeros(matT& mat, int N, int M = 0) 00111 { 00112 mat.resize(N, M == 0 ? N : M); 00113 mat = typename matT::value_type(0); 00114 } 00115 00119 template <class matT> 00120 void ones(matT& mat, int N, int M = 0) 00121 { 00122 mat.resize(N, M == 0 ? N : M); 00123 mat = typename matT::value_type(1); 00124 } 00125 00130 template <class matT> 00131 void eye(matT& mat, int N, int M = 0) 00132 { 00133 int nmin = (M == 0 ? N : (M < N ? M : N)); 00134 mat.resize(N, M == 0 ? N : M); 00135 mat = typename matT::value_type(0); 00136 typename matT::value_type one(1); 00137 for (int k = 0; k < nmin; ++k) 00138 mat(k, k) = one; 00139 } 00140 00148 template <class matT> 00149 void rand(matT &A, typename matT::value_type low = 0, 00150 typename matT::value_type high = 1) 00151 { 00152 int i, j, M = A.rows(), N = A.cols(); 00153 typename matT::value_type scale = high - low; 00154 for (j = 0; j < N; ++j) 00155 for (i = 0; i < M; ++i) 00156 A(i, j) = low + 00157 scale * double(std::rand()) / double(RAND_MAX); 00158 } 00159 00169 template <class matT> 00170 void int_rand(matT &A, typename matT::value_type low = 0, 00171 typename matT::value_type high = 1) 00172 { 00173 int i, j, M = A.rows(), N = A.cols(); 00174 double bins = high - low + 1; 00175 for (j = 0; j < N; ++j) 00176 for (i = 0; i < M; ++i) 00177 A(i, j) = low + 00178 typename matT::value_type( 00179 std::floor(bins * double(std::rand()) / 00180 double(RAND_MAX))); 00181 } 00182 00189 template <class matT> 00190 void from_diag(matT& mat, const matT& vect) 00191 { 00192 int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols()); 00193 mat = typename matT::value_type(0); 00194 if (vect.rows() != 1 && vect.cols() != 1) 00195 throw LaException("diag<matT>(const matT& vect, matT& mat)", 00196 "The argument 'vect' is not a vector: " 00197 "neither dimension is equal to one"); 00198 if (vect.rows() * vect.cols() != nmin) 00199 throw LaException("diag<matT>(const matT& vect, matT& mat)", 00200 "The size of the vector is unequal to the matrix diagonal"); 00201 if (vect.rows() == 1) 00202 for (int k = 0; k < nmin; ++k) 00203 mat(k, k) = vect(0, k); 00204 else 00205 for (int k = 0; k < nmin; ++k) 00206 mat(k, k) = vect(k, 0); 00207 } 00208 00210 00211 00217 template <class matT> 00218 matT zeros(int N, int M = 0) 00219 { 00220 matT mat; 00221 zeros(mat, N, M); 00222 return mat.shallow_assign(); 00223 } 00224 00228 template <class matT> 00229 matT ones(int N, int M = 0) 00230 { 00231 matT mat; 00232 ones(mat, N, M); 00233 return mat.shallow_assign(); 00234 } 00235 00239 template <class matT> 00240 matT eye(int N, int M = 0) 00241 { 00242 matT mat; 00243 eye(mat, N, M); 00244 return mat.shallow_assign(); 00245 } 00246 00255 template <class matT> 00256 matT rand(int N, int M, 00257 typename matT::value_type low = 0, 00258 typename matT::value_type high = 1) 00259 { 00260 matT mat(N, M); 00261 rand(mat, low, high); 00262 return mat.shallow_assign(); 00263 } 00264 00274 template <class matT> 00275 matT int_rand(int N, int M, 00276 typename matT::value_type low = 0, 00277 typename matT::value_type high = 1) 00278 { 00279 matT mat(N, M); 00280 int_rand(mat, low, high); 00281 return mat.shallow_assign(); 00282 } 00283 00289 template <class matT> 00290 matT from_diag(const matT& vect) 00291 { 00292 if (vect.rows() != 1 && vect.cols() != 1) 00293 throw LaException("diag<matT>(const matT& vect, matT& mat)", 00294 "The argument 'vect' is not a vector: " 00295 "neither dimension is equal to one"); 00296 int nmax(vect.rows() > vect.cols() ? vect.rows() : vect.cols()); 00297 matT mat(nmax, nmax); 00298 from_diag(mat, vect); 00299 return mat.shallow_assign(); 00300 } 00301 00318 template<class destT, class srcT> 00319 destT convert_mat(const srcT& src) 00320 { 00321 destT res(src.rows(), src.cols()); 00322 // optimize later; for now use the correct but slow implementation 00323 int i, j, M = src.rows(), N = src.cols(); 00324 for (j = 0; j < N; ++j) 00325 for (i = 0; i < M; ++i) 00326 res(i, j) = typename destT::value_type ( src(i, j) ); 00327 return res.shallow_assign(); 00328 } 00329 00333 template <class matT> 00334 matT linspace(typename matT::value_type start, 00335 typename matT::value_type end, 00336 int nr_points) 00337 { 00338 matT mat(nr_points, 1); 00339 typename matT::value_type stepsize = 00340 (end - start) / typename matT::value_type (nr_points - 1); 00341 for (int k = 0; k < nr_points; ++k) 00342 { 00343 mat(k, 0) = start; 00344 start += stepsize; 00345 } 00346 return mat.shallow_assign(); 00347 } 00348 00352 template <class matT> 00353 matT repmat(const matT& A, int M, int N) 00354 { 00355 int i, j, origM = A.rows(), origN = A.cols(); 00356 matT mat(origM * M, origN * N); 00357 for (j = 0; j < N; ++j) 00358 for (i = 0; i < M; ++i) 00359 mat(LaIndex(i * origM, (i + 1) * origM - 1), 00360 LaIndex(j * origN, (j + 1) * origN - 1)) 00361 .inject(A); 00362 00363 return mat.shallow_assign(); 00364 } 00366 00372 template <class matT> 00373 typename matT::value_type trace(const matT& mat) 00374 { 00375 int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols()); 00376 typename matT::value_type result(0); 00377 for (int k = 0; k < nmin; ++k) 00378 result += mat(k, k); 00379 return result; 00380 } 00381 00385 template <class matT> 00386 matT diag(const matT& mat) 00387 { 00388 int nmin = (mat.rows() < mat.cols() ? mat.rows() : mat.cols()); 00389 matT vect(nmin, 1); 00390 for (int k = 0; k < nmin; ++k) 00391 vect(k, 0) = mat(k, k); 00392 return vect.shallow_assign(); 00393 } 00394 00396 00397 } // namespace la 00398 00399 #endif // _LATMPL_H