90 if (
IsB ==
nullptr ) {
96 const PetscInt *loc_indices;
97 ISGetIndices(
IsA, &loc_indices);
99 a_used[ loc_indices[i] - ds->
begin()] =
true;
100 ISRestoreIndices(
IsA, &loc_indices);
104 PetscMalloc(
sizeof(PetscInt)*
loc_size_B, &b_indices);
106 if (! a_used[i]) b_indices[j++] = i + ds->
begin();
107 ISCreateGeneral(PETSC_COMM_WORLD,
loc_size_B, b_indices, PETSC_COPY_VALUES, &
IsB);
108 PetscFree(b_indices);
117 loc_size_A(other.loc_size_A), loc_size_B(other.loc_size_B), state(other.state),
118 Compl(other.Compl), ds_(other.ds_)
120 MatCopy(other.
A,
A, DIFFERENT_NONZERO_PATTERN);
121 MatCopy(other.
IA,
IA, DIFFERENT_NONZERO_PATTERN);
122 MatCopy(other.
IAB,
IAB, DIFFERENT_NONZERO_PATTERN);
171 PetscErrorCode ierr = 0;
173 MatStructure mat_subset_pattern;
174 PetscScalar *sol_array;
176 mat_reuse=MAT_REUSE_MATRIX;
177 mat_subset_pattern=SUBSET_NONZERO_PATTERN;
179 mat_reuse=MAT_INITIAL_MATRIX;
180 mat_subset_pattern=DIFFERENT_NONZERO_PATTERN;
185 VecCreateMPI(PETSC_COMM_WORLD,
loc_size_A, PETSC_DETERMINE, &(
RHS1));
186 VecCreateMPI(PETSC_COMM_WORLD,
loc_size_B, PETSC_DETERMINE, &(
RHS2));
191 VecCreateMPI(PETSC_COMM_WORLD,
loc_size_A, PETSC_DETERMINE, &(
Sol1));
192 VecCreateMPI(PETSC_COMM_WORLD,
loc_size_B, PETSC_DETERMINE, &(
Sol2));
196 VecGetArray(
Sol2, &sol_array );
198 VecRestoreArray(
Sol2, &sol_array );
226 ierr+=MatMatMult(
IA,
B, mat_reuse, 1.0 ,&(
IAB));
229 ierr+=MatMatMult(
Bt,
IAB, mat_reuse, 1.9 ,&(
xA));
247 ASSERT( ierr == 0 ).error(
"PETSC Error during calculation of Schur complement.\n");
262 VecScatterBegin(
rhs1sc,
rhs_,
RHS1, INSERT_VALUES, SCATTER_FORWARD);
263 VecScatterEnd(
rhs1sc,
rhs_,
RHS1, INSERT_VALUES, SCATTER_FORWARD);
264 VecScatterBegin(
rhs2sc,
rhs_,
RHS2, INSERT_VALUES, SCATTER_FORWARD);
265 VecScatterEnd(
rhs2sc,
rhs_,
RHS2, INSERT_VALUES, SCATTER_FORWARD);
296 if (
ds_ != NULL)
delete ds_;
304 PetscInt ncols, pos_start, pos_start_IA;
306 MatReuse mat_reuse=MAT_REUSE_MATRIX;
310 MatDuplicate(
A, MAT_DO_NOT_COPY_VALUES, &
IA);
313 MatGetOwnershipRange(
A,&pos_start,PETSC_NULL);
314 MatGetOwnershipRange(
IA,&pos_start_IA,PETSC_NULL);
317 const PetscInt *
cols;
318 const PetscScalar *vals;
322 unsigned int mat_block=1;
323 for(
unsigned int loc_row=0; loc_row < processed_rows.size(); loc_row++) {
324 if (processed_rows[loc_row] != 0)
continue;
326 PetscInt min=std::numeric_limits<int>::max(), max=-1, size_submat;
329 MatGetRow(
A, loc_row + pos_start, &ncols, &
cols, PETSC_NULL);
330 for (PetscInt i=0; i<ncols; i++) {
342 size_submat = max - min + 1;
343 ASSERT(ncols-b_vals == size_submat).error(
"Submatrix cannot contains empty values.\n");
345 MatRestoreRow(
A, loc_row + pos_start, &ncols, &
cols, PETSC_NULL);
346 arma::mat submat2(size_submat, size_submat);
348 for (PetscInt i=0; i<size_submat; i++) {
349 processed_rows[ loc_row + i ] = mat_block;
350 submat_rows.push_back( i + loc_row + pos_start_IA );
351 MatGetRow(
A, i + loc_row + pos_start, &ncols, &
cols, &vals);
352 for (PetscInt j=0; j<ncols; j++) {
354 submat2( i,
cols[j] - loc_row - pos_start ) = vals[j];
357 MatRestoreRow(
A, i + loc_row + pos_start, &ncols, &
cols, &vals);
362 const PetscInt*
rows = &submat_rows[0];
363 MatSetValues(
IA, submat_rows.size(),
rows, submat_rows.size(),
rows, invmat.memptr(), INSERT_VALUES);
368 MatAssemblyBegin(
IA, MAT_FINAL_ASSEMBLY);
369 MatAssemblyEnd(
IA, MAT_FINAL_ASSEMBLY);
383 return std::numeric_limits<double>::infinity();
422 START_TIMER(
"SchurComplemet::resolve without form schur");
448 if (
A != NULL )
chkerr(MatDestroy(&
A));
449 if (
B != NULL )
chkerr(MatDestroy(&
B));
451 if (
C != NULL )
chkerr(MatDestroy(&
C));
468 if (
ds_ != NULL)
delete ds_;
#define ASSERT_PTR(ptr)
Definition of assert macro checking non-null pointer (PTR) only for debug mode.
unsigned int begin(int proc) const
get starting local index
unsigned int lsize(int proc) const
get local size
double compute_residual() override
Mat matrix_
Petsc matrix of the problem.
LinSys::SolveInfo solve() override
void set_from_input(const Input::Record in_rec) override
Vec rhs_
PETSc vector constructed with vx array.
const Mat * get_matrix() override
double get_solution_precision() override
const Vec * get_rhs() override
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
const Vec & get_solution()
void set_solution(Vec sol_vec)
friend class SchurComplement
void set_matrix_changed()
bool is_negative_definite()
Vec solution_
PETSc vector constructed with vb array.
bool matrix_changed_
true if the matrix was changed since the last solve
bool rhs_changed_
true if the right hand side was changed since the last solve
Distribution * make_complement_distribution()
get distribution of complement object if complement is defined
void set_from_input(const Input::Record in_rec) override
void set_tolerances(double r_tol, double a_tol, double d_tol, unsigned int max_it) override
double compute_residual() override
void set_complement(LinSys_PETSC *ls)
Set complement LinSys object.
LinSys::SolveInfo solve() override
double get_solution_precision() override
get precision of solving
void create_inversion_matrix()
create IA matrix
Support classes for parallel programing.
Wrappers for linear systems based on MPIAIJ and MATIS format.
#define DebugOut()
Macro defining 'debug' record of log.
ArmaMat< double, N, M > mat
Assembly explicit Schur complement for the given linear system. Provides method for resolution of the...
#define START_TIMER(tag)
Starts a timer with specified tag.
void chkerr(unsigned int ierr)
Replacement of new/delete operator in the spirit of xmalloc.