bareiss_test.cpp

/******************************************************************************
* MODULE     : bareiss_test.cpp
* DESCRIPTION: Test matrices
* COPYRIGHT  : (C) 2008  Bernard Mourrain
*******************************************************************************
* This software falls under the GNU general public license and comes WITHOUT
* ANY WARRANTY WHATSOEVER. See the file $TEXMACS_PATH/LICENSE for more details.
* If you don't have this file, write to the Free Software Foundation, Inc.,
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
******************************************************************************/

#include <numerix/integer.hpp>
#include <numerix/rational.hpp>
#include <algebramix/matrix_bareiss.hpp>
using namespace mmx;

void
bareiss_test () {
  int t[]={1,2,3,4,5,6,7,8,9,10,11,12};
  integer* o1= mmx_new<integer> (9);
  for (nat i = 0; i < 9; i++) o1[i] = integer (t[i]);
  matrix<integer> A1(o1,3,3,format<integer>());
  mmout << "A1 = " << A1 << "\n";
  mmout << "rk = " << bareiss_rank(A1) << "\n";
  mmout << "ker= " << bareiss_kernel(A1) << "\n";
  mmout << "chk= " << A1*bareiss_kernel(A1) << "\n";
  mmout << "cok= " << bareiss_cokernel(A1) << "\n";
  mmout << "chk= " << bareiss_cokernel(A1)*A1 << "\n";
  mmout << "img= " << bareiss_image(A1) << "\n";
  mmout << "coi= " << bareiss_coimage(A1) << "\n";
  bareiss_pivoting (A1);
  mmout << "A1t= " << A1 << "\n";

  integer* o2= mmx_new<integer> (12);
  for (nat i = 0; i < 12; i++) o2[i] = integer (t[i]);
  matrix<integer> A2(o2,4,3,format<integer>());
  mmout << "A2 = " << A2 << "\n";
  //  unsigned r;
  mmout << "rk = " << bareiss_rank(A2) << "\n";
  mmout << "ker= " << bareiss_kernel(A2) << "\n";
  mmout << "chk= " << A2*bareiss_kernel(A2) << "\n";
  mmout << "cok= " << bareiss_cokernel(A2) << "\n";
  mmout << "chk= " << bareiss_cokernel(A2)*A2 << "\n";
  mmout << "img= " << bareiss_image(A2) << "\n";
  mmout << "coi= " << bareiss_coimage(A2) << "\n";
  bareiss_pivoting (A2);
  mmout << "A2t= " << A2 << "\n";

  integer* o3= mmx_new<integer> (12);
  for (nat i = 0; i < 12; i++) o3[i] = integer (t[i]);
  matrix<integer> A3(o3,3,4,format<integer>());
  mmout << "A3 = " << A3 << "\n";
  //  unsigned r;
  mmout << "rk = " << bareiss_rank(A3) << "\n";
  mmout << "ker= " << bareiss_kernel(A3) << "\n";
  mmout << "chk= " << A3*bareiss_kernel(A3) << "\n";
  mmout << "cok= " << bareiss_cokernel(A3) << "\n";
  mmout << "chk= " << bareiss_cokernel(A3)*A3 << "\n";
  mmout << "img= " << bareiss_image(A3) << "\n";
  mmout << "coi= " << bareiss_coimage(A3) << "\n";
  bareiss_pivoting (A3);
  mmout << "A3t= " << A3 << "\n";

  integer* m= mmx_new<integer> (400);
  integer* n= m;
  for (nat i = 0; i < 400; i++) {
    *n = integer (i+1);
    n++;
  }
  matrix<integer> U (m,20,20,format<integer>()) ;
  mmout << "U = \n"<< U << "\n";
  bareiss_pivoting (U);
  mmout << "U = \n"<< U << "\n";

  rational* a= mmx_new<rational> (400);
  rational* b= a;
  for (nat i = 0; i < 400; i++) {
    *b = rational (i+1);
    b++;
  }
  matrix<rational> Ar (a,20,20,format<rational>());
  mmout << "det= " << bareiss_det (Ar)<<"\n";

//   matrix<rational> B (rational(0), 5, 3);
//   B(0,0)= 1; B(0,1)= 2;
//   B(1,0)= 7; B(1,1)= 3;
//   B(2,0)= 6; B(2,1)= 5;
//   B(4,2)= 11;
//   mmout << "B= " << B << "\n";
//   mmout << "ker(B)= " << kernel (B) << "\n";
//   mmout << "ker(B)B= " << kernel (B) * B << "\n";
}

int
main() {
  try {
    bareiss_test ();
  }
  catch (exception& e) {
    mmout << "ERROR in pivot_test: " << e << "\n";
    return -1;
  }
  return 0;
}