dolfin team mailing list archive

["Alessio Quaglino" <alessio@xxxxxx>]

I've discovered that the unusual slowness of umfpack does not depend on
windows (and thus I finally managed to compile and link everything, ATLAS
BLAS included) but on the following problem that I could replicate also in
linux (ubuntu).

When I use the standard uBlasLUSolver class everything works at the
expected speed, while when I just copy and paste the same code in my .cpp
solver it becomes about 50-100 times slower. However, I cannot use the
dolfin uBlasLUSolver, since I need to save the numeric objects and use
them several time, and I need to modify some UMFPACK options as well.

My class looks like this:

//-----------------------------------------------------------------------------
MyLUSolver::MyLUSolver(char* number) /*: uBlasLUSolver()*/
{
  fact = 0;
  tag = number;
}
//-----------------------------------------------------------------------------
void MyLUSolver::factorize(dolfin::uint flag)
{
  fact = flag;
}
//-----------------------------------------------------------------------------
uint MyLUSolver::solve(const uBlasSparseMatrix& A, uBlasVector& x, const
uBlasVector& b)
{
  // Check dimensions and get number of non-zeroes
  const uint M  = A.size(0);
  const uint N  = A.size(1);
  const uint nz = A.nnz();

  dolfin_assert(M == A.size(1));
  dolfin_assert(nz > 0);

  x.init(N);
  message("Solving linear system of size %d x %d (UMFPACK LU solver).", M,
N);

  // Make sure matrix assembly is complete
  (const_cast< uBlasSparseMatrix& >(A)).complete_index1_data();

  int status;
  long int* inull = (long int *) NULL;
  void *Symbolic, *Numeric;
  const std::size_t* Ap = &(A.index1_data() [0]);
  const std::size_t* Ai = &(A.index2_data() [0]);
  const double* Ax = &(A.value_data() [0]);
  double* xx = &(x.data() [0]);
  const double* bb = &(b.data() [0]);

  // Set controls
  double Control [UMFPACK_CONTROL];
  Control [UMFPACK_PRL] = 3;
  double Info [UMFPACK_INFO];

  // Solve for transpose since UMFPACK expects compressed column format
long int* Rp = new long int[M+1];
  long int* Ri = new long int[nz];
  double* Rx   = new double[nz];

  status = umfpack_dl_transpose(M, M, (const long int*) Ap, (const long
int*) Ai, Ax, inull, inull, Rp, Ri, Rx);
  if (status < 0)
	  std::cerr << "umfpack_dl_transpose failed. Error code: " << status <<
std::endl;

  // Init files
  char fs[20];
  char fn[20];
  strcpy_s(fs, tag);
  strcat_s(fs, "_A2LUsym.bin");
  strcpy_s(fn, tag);
  strcat_s(fn, "_A2LUnum.bin");

  if( fact == 0 )
  {
    // Symbolic factorization
    status = umfpack_dl_symbolic(M, M, (const long int*) Rp, (const long
int*) Ri, Rx, &Symbolic, Control, Info);
    //umfpack_dl_report_info(Control, Info);

	if (status < 0)
        std::cerr << "umfpack_dl_symbolic failed. Error code: " << status
<< std::endl;

	// Save symbolic object
    status = umfpack_dl_save_symbolic(Symbolic, fs);
    if (status < 0)
        std::cerr << "umfpack_dl_save_symbolic failed. Error code: " <<
status << std::endl;

	// Clean memory
    umfpack_dl_free_symbolic(&Symbolic);
  }

  else if( fact == 1 )
  {
	// Load symbolic object
    status = umfpack_dl_load_symbolic(&Symbolic, fs);
    if (status < 0)
        std::cerr << "umfpack_dl_load_symbolic failed. Error code: " <<
status << std::endl;

	// Numeric factorization
    umfpack_dl_numeric( (const long int*) Rp, (const long int*) Ri, Rx,
Symbolic, &Numeric, Control, Info);
	//umfpack_dl_report_info(Control, Info);
    umfpack_dl_free_symbolic(&Symbolic);

	// Save numeric object
    status = umfpack_dl_save_numeric(Numeric, fn);
    if (status < 0)
        std::cerr << "umfpack_dl_save_numeric failed. Error code: " <<
status << std::endl;

	// Clean memory
    umfpack_dl_free_numeric(&Numeric);
  }

  else
  {
    // Load numeric object
    status = umfpack_dl_load_numeric(&Numeric, fn);
    if (status < 0)
        std::cerr << "umfpack_di_load_numeric failed. Error code: " <<
status << std::endl;

	// Solve linear system
    umfpack_dl_solve(UMFPACK_A, (const long int*) Rp, (const long int*)
Ri, Rx, xx, bb, Numeric, Control, Info);

	// Claen memory
    umfpack_dl_free_numeric(&Numeric);
  }

  delete [] Rp;
  delete [] Ri;
  delete [] Rx;
  return 1;
}
//-----------------------------------------------------------------------------

Alessio