// $Id$
// Version: $Name$
//
-// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010, 2011 by the deal.II authors
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2009, 2010, 2011, 2012 by the deal.II authors
//
// This file is subject to QPL and may not be distributed
// without copyright and license information. Please refer
"re-configuring the library and re-building it!"));
}
-
+
extern "C"
void ma27bd_ (const unsigned int *,
const unsigned int *,
"the respective files in the right place, "
"re-configuring the library and re-building it!"));
}
-
+
extern "C" void ma27x2_ (unsigned int *)
{
"the respective files in the right place, "
"re-configuring the library and re-building it!"));
}
-
-
+
+
extern "C" void ma27x3_ (const unsigned int *)
{
AssertThrow (false,
"the respective files in the right place, "
"re-configuring the library and re-building it!"));
}
-
+
extern "C"
void ma47ad_ (const unsigned int *,
"re-configuring the library and re-building it!"));
}
-
+
extern "C"
void ma47bd_ (const unsigned int *,
const unsigned int *,
"re-configuring the library and re-building it!"));
}
-
+
extern "C"
void ma47cd_ (const unsigned int *,
const double *,
-namespace
+namespace
{
/**
* Output an error message and terminate the program.
Threads::ThreadMutex SparseDirectMA27::static_synchronisation_lock;
-struct SparseDirectMA27::DetachedModeData
+struct SparseDirectMA27::DetachedModeData
{
/**
* Mutex to assure that only one
* through the pipe.
*/
Threads::ThreadMutex mutex;
-
+
/**
* File handles for the pipe
* between server (computing
count += ret;
};
-
+
std::fflush (NULL);
}
-
+
/**
* Get a message from the client
* program. Obey all the rules
reinterpret_cast<char *> (t) + count,
sizeof(T) * N - count);
while ((ret<0) && (errno==EINTR));
-
+
if (ret < 0)
die ("error on client side in 'get'", ret, errno, child_pid);
-SparseDirectMA27::~SparseDirectMA27()
+SparseDirectMA27::~SparseDirectMA27()
{
if (detached_mode)
if (detached_mode_data != 0)
void
-SparseDirectMA27::set_detached_mode ()
+SparseDirectMA27::set_detached_mode ()
{
Assert (initialize_called == false,
ExcInitializeAlreadyCalled());
// Assign the return value to a
// variable to avoid compiler
// warnings
-//TODO:[WB] Use t to trace errors?
+//TODO:[WB] Use t to trace errors?
int t = pipe(detached_mode_data->server_client_pipe);
- (void)t;
+ (void)t;
// fflush(NULL) is said to be a
// good idea before fork()
const pid_t parent_pid = getpid();
detached_mode_data->put (&parent_pid, 1, "parent_pid");
};
-
-
+
+
// suppress error output if
// requested
if (suppress_output)
const unsigned int LP = 0;
call_ma27x3 (&LP);
};
-
+
sparsity_pattern = &sp;
-
+
const unsigned int
n_rows = sparsity_pattern->n_rows();
const std::size_t * const
++col)
if (row <= *col)
++n_nonzero_elements;
-
+
// fill the row numbers and column
// numbers arrays from the sparsity
if (row <= *col)
{
Assert (global_index < n_nonzero_elements, ExcInternalError());
-
+
row_numbers[global_index] = row+1;
column_numbers[global_index] = *col+1;
++global_index;
};
Assert (global_index == n_nonzero_elements, ExcInternalError());
-
+
// initialize scratch arrays and
// variables
LIW = static_cast<unsigned int>((2*n_nonzero_elements + 3*n_rows + 1) *
// allowed to allocate more memory
// no more
bool call_succeeded = true;
- do
+ do
{
call_ma27ad (&n_rows, &n_nonzero_elements,
&row_numbers[0], &column_numbers[0],
// increase allowed: exit loop
if (call_succeeded || (LIW_increase_factor_1 <= 1))
break;
-
+
// otherwise: increase LIW and retry
LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_1);
IW.resize (LIW);
// are the same
Assert (sparsity_pattern == &matrix.get_sparsity_pattern(),
ExcDifferentSparsityPatterns());
-
-
+
+
// set LA and fill the A array of
// values
LA = std::max (static_cast<int>(NRLNEC * LA_factor),
LIW = static_cast<unsigned int>(NIRNEC * LIW_factor_2);
IW.resize (LIW);
};
-
+
const unsigned int n_rows = matrix.get_sparsity_pattern().n_rows();
-
+
// loop until memory requirements
// are satisfied or we are not
// allowed to allocate more memory
// no more
bool call_succeeded = true;
- do
+ do
{
call_ma27bd (&n_rows, &n_nonzero_elements,
&row_numbers[0], &column_numbers[0],
break;
};
-
+
// otherwise: increase LIW or
// LA if that is allowed and
{
if (LIW_increase_factor_2 <= 1)
goto exit_loop;
-
+
LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_2);
IW.resize (LIW);
break;
// an exception on the
// allocation.
std::cout << "<*>" << std::flush;
-
+
LA = static_cast<unsigned int>(LA * LA_increase_factor);
if (true)
{
std::vector<double> tmp;
A.swap (tmp);
};
-
+
A.resize (LA);
fill_A (matrix);
-
+
break;
};
-
+
// ups, other return
// value, don't know
// what to do here
SparseDirectMA27::solve (Vector<double> &rhs_and_solution) const
{
Assert (factorize_called == true, ExcFactorizeNotCalled());
-
+
const unsigned int n_rows = rhs_and_solution.size();
call_ma27cd (&n_rows, &A[0], &LA,
&IW[0], &LIW, &MAXFRT,
// doubles
Vector<double> tmp (rhs_and_solution.size());
tmp = rhs_and_solution;
-
+
const unsigned int n_rows = rhs_and_solution.size();
call_ma27cd (&n_rows, &A[0], &LA,
&IW[0], &LIW, &MAXFRT,
if (detached_mode)
return non_static_synchronisation_lock;
else
- return static_synchronisation_lock;
+ return static_synchronisation_lock;
}
Assert (n_nonzero_elements <= A.size(), ExcInternalError());
const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
-
+
const unsigned int n_rows = sparsity_pattern.n_rows();
const std::size_t *rowstart_indices = sparsity_pattern.get_rowstart_indices();
const unsigned int *col_nums = sparsity_pattern.get_column_numbers();
if (row <= *col)
{
Assert (global_index < n_nonzero_elements, ExcInternalError());
-
+
A[global_index] = matrix(row,*col);
++global_index;
(std::fabs(matrix(row,*col) - matrix(*col,row))
<= 1e-15 * std::fabs (matrix(row,*col))),
ExcMatrixNotSymmetric());
-
- Assert (global_index == n_nonzero_elements, ExcInternalError());
+
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
}
-
+
void SparseDirectMA27::call_ma27ad (const unsigned int *N,
const unsigned int *NZ,
const unsigned int *IRN,
// function index, then array
// sizes, then arrays
detached_mode_data->put ("1", 1, "ACTION 1");
-
+
detached_mode_data->put (N, 1, "N");
detached_mode_data->put (NZ, 1, "NZ");
detached_mode_data->put (IRN, *NZ, "IRN");
// be in used on this side
Assert (this->IKEEP.size() == 0, ExcInternalError());
Assert (this->IW1.size() == 0, ExcInternalError());
-
+
// next get back what we need
// to know
detached_mode_data->get (IFLAG, 1, "IFLAG");
// except for A and LA
Threads::ThreadMutex::ScopedLock lock (detached_mode_data->mutex);
detached_mode_data->put ("2", 1, "ACTION 2");
-
+
detached_mode_data->put (LA, 1, "LA");
detached_mode_data->put (A, *LA, "A");
};
}
-
+
// then start with work
matrix = &m;
const SparsityPattern &sparsity_pattern = matrix->get_sparsity_pattern();
-
+
const unsigned int
n_rows = sparsity_pattern.n_rows();
const std::size_t * const
// required by the docs of MA47
if ((row <= *col) && (m(row,*col) != 0))
++n_nonzero_elements;
-
+
// fill the row numbers and column
// numbers arrays from the sparsity
if ((row <= *col) && (m(row,*col) != 0))
{
Assert (global_index < n_nonzero_elements, ExcInternalError());
-
+
row_numbers[global_index] = row+1;
column_numbers[global_index] = *col+1;
++global_index;
};
Assert (global_index == n_nonzero_elements, ExcInternalError());
-
+
// initialize scratch arrays and
// variables
LIW = static_cast<unsigned int>((2*n_nonzero_elements + 5*n_rows + 4) *
// increase allowed: exit loop
if (call_succeeded || (LIW_increase_factor_1 <= 1))
break;
-
+
// otherwise: increase LIW and retry
LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_1);
- IW.resize (LIW);
+ IW.resize (LIW);
}
while (true);
{
Assert (factorize_called == false,
ExcCantFactorizeAgain());
-
+
// if necessary, initialize process
if (initialize_called == false)
initialize (m);
// make sure the matrices
// are the same
Assert (matrix == &m, ExcDifferentMatrices());
-
-
+
+
// set LA and fill the A array of
// values
LA = std::max (static_cast<int>(INFO[5] * LA_factor),
static_cast<int>(n_nonzero_elements));
A.resize (LA);
fill_A (m);
-
+
// if necessary extend IW
if (LIW < INFO[6] * LIW_factor_2)
{
// output info flags
bool call_succeeded;
- do
+ do
{
call_ma47bd (&n_rows, &n_nonzero_elements, &column_numbers[0],
&A[0], &LA,
{
if (LIW_increase_factor_2 <= 1)
goto exit_loop;
-
+
LIW = static_cast<unsigned int>(LIW * LIW_increase_factor_2);
IW.resize (LIW);
break;
// an exception on the
// allocation.
std::cout << "<*>" << std::flush;
-
+
LA = static_cast<unsigned int>(LA * LA_increase_factor);
if (true)
{
std::vector<double> tmp;
A.swap (tmp);
};
-
+
A.resize (LA);
fill_A (m);
-
+
break;
};
-
+
// ups, other return
// value, don't know
// what to do here
SparseDirectMA47::solve (Vector<double> &rhs_and_solution)
{
Assert (factorize_called == true, ExcFactorizeNotCalled());
-
+
const unsigned int n_rows = rhs_and_solution.size();
call_ma47cd (&n_rows, &A[0], &LA,
&IW[0], &LIW,
Assert (n_nonzero_elements <= A.size(), ExcInternalError());
const SparsityPattern &sparsity_pattern = matrix.get_sparsity_pattern ();
-
+
const unsigned int n_rows = sparsity_pattern.n_rows();
const std::size_t *rowstart_indices = sparsity_pattern.get_rowstart_indices();
const unsigned int *col_nums = sparsity_pattern.get_column_numbers();
if ((row <= *col) && (matrix(row,*col) != 0))
{
Assert (global_index < n_nonzero_elements, ExcInternalError());
-
+
A[global_index] = matrix(row,*col);
++global_index;
||
(matrix(row,*col) == matrix(*col,row)),
ExcMatrixNotSymmetric());
-
- Assert (global_index == n_nonzero_elements, ExcInternalError());
+
+ Assert (global_index == n_nonzero_elements, ExcInternalError());
}
const unsigned int *ICNTL, //length 7
int *INFO) //length 24
{
- double RINFO[4];
+ double RINFO[4];
HSL::MA47::ma47ad_(n_rows, n_nonzero_elements,
row_numbers, column_numbers,
IW, LIW, KEEP,
unsigned int *IW1, //2*n_rows+2
int *INFO) //length 24
{
- double RINFO[4];
+ double RINFO[4];
HSL::MA47::ma47bd_(n_rows, n_nonzero_elements, column_numbers,
A, LA,
IW, LIW, KEEP, CNTL, ICNTL,
std::vector<double> W(*n_rows);
HSL::MA47::ma47cd_(n_rows, A, LA,
IW, LIW, &W[0],
- rhs_and_solution, IW1, ICNTL);
+ rhs_and_solution, IW1, ICNTL);
}
umfpack_dl_free_numeric (&numeric_decomposition);
numeric_decomposition = 0;
}
-
+
{
std::vector<long int> tmp;
tmp.swap (Ap);
std::vector<long int> tmp;
tmp.swap (Ai);
}
-
+
{
std::vector<double> tmp;
tmp.swap (Ax);
factorize (const Matrix &matrix)
{
Assert (matrix.m() == matrix.n(), ExcNotQuadratic())
-
+
clear ();
const unsigned int N = matrix.m();
Ap[row] = Ap[row-1] + matrix.get_row_length(row-1);
Assert (static_cast<unsigned int>(Ap.back()) == Ai.size(),
ExcInternalError());
-
+
// then copy over matrix
// elements. note that for sparse
// matrices, iterators are sorted
// row points to the first entry
// not yet written to
std::vector<long int> row_pointers = Ap;
-
+
for (typename Matrix::const_iterator p=matrix.begin();
p!=matrix.end(); ++p)
{
// matrices, so ship this task out
// to a different function
sort_arrays (matrix);
-
+
int status;
status = umfpack_dl_symbolic (N, N,
&Ap[0], &Ai[0], &Ax[0],
&control[0], 0);
AssertThrow (status == UMFPACK_OK,
ExcUMFPACKError("umfpack_dl_symbolic", status));
-
+
status = umfpack_dl_numeric (&Ap[0], &Ai[0], &Ax[0],
symbolic_decomposition,
&numeric_decomposition,
Assert (Ap.size() != 0, ExcNotInitialized());
Assert (Ai.size() != 0, ExcNotInitialized());
Assert (Ai.size() == Ax.size(), ExcNotInitialized());
-
+
Vector<double> rhs (rhs_and_solution.size());
rhs = rhs_and_solution;
-
+
// solve the system. note that since
// UMFPACK wants compressed column storage
// instead of the compressed row storage
template <class Matrix>
void SparseDirectUMFPACK::factorize (const Matrix &)
{
- AssertThrow(false, ExcNeedsUMFPACK());
+ AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but you called ./configure without the necessary --with-umfpack switch."));
}
void
SparseDirectUMFPACK::solve (Vector<double> &) const
{
- AssertThrow(false, ExcNeedsUMFPACK());
+ AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but you called ./configure without the necessary --with-umfpack switch."));
}
SparseDirectUMFPACK::solve (const Matrix &,
Vector<double> &)
{
- AssertThrow(false, ExcNeedsUMFPACK());
+ AssertThrow(false, ExcMessage("To call this function you need UMFPACK, but you called ./configure without the necessary --with-umfpack switch."));
}
}
#ifdef DEAL_II_USE_MUMPS
-SparseDirectMUMPS::SparseDirectMUMPS ()
+SparseDirectMUMPS::SparseDirectMUMPS ()
:
initialize_called (false)
{}
-SparseDirectMUMPS::~SparseDirectMUMPS ()
+SparseDirectMUMPS::~SparseDirectMUMPS ()
{}
template <class Matrix>
-void SparseDirectMUMPS::initialize (const Matrix& matrix,
- const Vector<double> & vector)
+void SparseDirectMUMPS::initialize (const Matrix& matrix,
+ const Vector<double> & vector)
{
// Check we haven't been here before:
Assert (initialize_called == false, ExcInitializeAlreadyCalled());
// Use MPI_COMM_WORLD as communicator
id.comm_fortran = -987654;
dmumps_c (&id);
-
+
// Hand over matrix and right-hand side
- if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+ if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
{
// Objects denoting a MUMPS data structure:
// number of nonzero elements in matrix
nz = matrix.n_actually_nonzero_elements ();
- // representation of the matrix
+ // representation of the matrix
a = new double[nz];
// matrix indices pointing to the row and
// is the matrix element (irn[k], jcn[k])
irn = new unsigned int[nz];
jcn = new unsigned int[nz];
-
+
unsigned int index = 0;
-
- for (SparseMatrix<double>::const_iterator ptr = matrix.begin ();
+
+ for (SparseMatrix<double>::const_iterator ptr = matrix.begin ();
ptr != matrix.end (); ++ptr)
- if (std::abs (ptr->value ()) > 0.0)
+ if (std::abs (ptr->value ()) > 0.0)
{
a[index] = ptr->value ();
irn[index] = ptr->row () + 1;
jcn[index] = ptr->column () + 1;
++index;
}
-
+
rhs = new double[n];
-
+
for (unsigned int i = 0; i < n; ++i)
rhs[i] = vector (i);
-
+
id.n = n;
id.nz = nz;
id.irn = irn;
id.a = a;
id.rhs = rhs;
}
-
- // No outputs
+
+ // No outputs
id.icntl[0] = -1;
id.icntl[1] = -1;
id.icntl[2] = -1;
initialize_called = true;
}
-void SparseDirectMUMPS::solve (Vector<double>& vector)
+void SparseDirectMUMPS::solve (Vector<double>& vector)
{
// Check that the solver has been initialized
// by the routine above:
dmumps_c (&id);
// Copy solution into the given vector
- if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
+ if (Utilities::MPI::this_mpi_process (MPI_COMM_WORLD) == 0)
{
for (unsigned int i=0; i<n; ++i)
vector(i) = rhs[i];
-
+
delete[] a;
delete[] irn;
delete[] jcn;
}
#endif // DEAL_II_USE_MUMPS
-// explicit instantiations for SparseMatrixMA27
+// explicit instantiations for SparseMatrixMA27
template
void SparseDirectMA27::factorize (const SparseMatrix<double> &matrix);
// explicit instantiations for SparseDirectMUMPS
#ifdef DEAL_II_USE_MUMPS
template <class Matrix>
-void SparseDirectMUMPS::initialize (const Matrix& matrix,
+void SparseDirectMUMPS::initialize (const Matrix& matrix,
const Vector<double> & vector);
-void SparseDirectMUMPS::solve (Vector<double>& vector);
+void SparseDirectMUMPS::solve (Vector<double>& vector);
#endif
DEAL_II_NAMESPACE_CLOSE