Lapack++
latmpl.h
Go to the documentation of this file.
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