"Fossies" - the Fresh Open Source Software Archive

Member "gama-2.05/tests/matvec/matvec_test_004.cpp" (10 May 2019, 3209 Bytes) of package /linux/privat/gama-2.05.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file.

    1 /* matvec_test_004.cpp
    2    Copyright (C) 2000, 2012  Ales Cepek <cepek@gnu.org>
    3 
    4    This library is free software; you can redistribute it and/or
    5    modify it under the terms of the GNU Library General Public
    6    License as published by the Free Software Foundation; either
    7    version 2 of the License, or (at your option) any later version.
    8 
    9    This library is distributed in the hope that it will be useful,
   10    but WITHOUT ANY WARRANTY; without even the implied warranty of
   11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
   12    Library General Public License for more details.
   13 
   14    You should have received a copy of the GNU Library General Public
   15    License along with this library (see COPYING.LIB); if not, write to the
   16    Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
   17 */
   18 
   19 #include <matvec/matvec.h>
   20 #include <matvec/pinv.h>
   21 #include <iostream>
   22 #include <iomanip>
   23 
   24 using namespace GNU_gama;
   25 
   26 template <class Float, class Exc> void
   27 InitMat_CACM_Alg_52(int n, Mat<Float, Exc>& A, Mat<Float, Exc>& invA)
   28 {
   29   /* CACM Algorithm 52
   30    * A set of test matrices by J.R.Herndon, corrected by P. Naur
   31    * -------------------------------------------------------------------
   32    *
   33    * This procedure places in A an n x n matrix whose inverse and
   34    * eigenvalues are known. The n-th row and the n-th column of the
   35    * inverse are the set: 1, 2, 3, ..., n. The matrix formed by
   36    * deleting the n-th row and n-th column of the inverse is the
   37    * identity matrix of order n-1.
   38    *
   39    */
   40 
   41   {
   42     A.reset(n,n);
   43 
   44     Float c = n * (n+1) * (n+n-5)/6;
   45     Float d = 1 / c;
   46 
   47     A(n,n) = -d;
   48     for (int i=1; i<=n-1; i++)
   49       {
   50         A(i,n) = A(n,i) = d * i;
   51         A(i,i) = d * (c - i*i);
   52 
   53         for (int j=1; j<=i-1; j++)
   54           A(i,j) = A(j,i) = -d*i*j;
   55       }
   56   }
   57 
   58   {
   59     invA.reset(n,n);
   60     invA.set_identity();
   61     for (int i=1; i<=n; i++) invA(i,n) = invA(n,i) = i;
   62   }
   63 }
   64 
   65 
   66 template <class Float, class Exc>
   67 Float MaxEl(const GNU_gama::MatVecBase<Float, Exc>& M)
   68 {
   69   Float m = 0, x;
   70   typedef typename GNU_gama::MatVecBase<Float, Exc>::const_iterator cIter;
   71   for (cIter i=M.begin(); i!=M.end(); ++i)
   72     {
   73       x = *i;  if (x < 0) x = -x;
   74       if (x > m) m = x;
   75     }
   76   return m;
   77 }
   78 
   79 int main()
   80 {
   81   using namespace std;
   82   using namespace GNU_gama;
   83 
   84   cout << "\n   CACM Alg. 52  ......   test_004  matvec "
   85        << GNU_gama::matvec_version() << "\n"
   86        << "------------------------------------------------------\n\n";
   87 
   88   const int MaxDim=70;
   89 
   90   for (int N=3; N<=MaxDim; N++)
   91     {
   92       cout << N << "x" << N;
   93 
   94       Mat<> A(N, N), B(N,N), C;
   95       InitMat_CACM_Alg_52(N, A, B);
   96 
   97       cout.precision(2);
   98       try {
   99         Mat<> I(N,N);
  100         I.set_identity();
  101 
  102         cout << "\tpinv() ";
  103         C = pinv(A);
  104         cout << MaxEl(B-C);
  105 
  106         cout << "\tinvert() ";
  107         C = A;
  108         C.invert();
  109         cout << MaxEl(B-C);
  110       }
  111       catch (Exception::matvec e) {
  112         cout << " err = " << int(e.error) << " : " << e.description;
  113       }
  114       catch (...) {
  115         cout << "?????????????????????????????????\n";
  116       }
  117       cout << endl;
  118     }
  119 
  120   cout <<  "\n------------------------------------------------------\n\n";
  121  }