47 vec_ds(vec_lsize, PETSC_COMM_WORLD),
49 positive_definite(false),
53 v_rhs=(
double *) xmalloc(
sizeof(
double) * (this->
vec_lsize() + 1) );
54 VecCreateMPIWithArray(PETSC_COMM_WORLD,1, this->
vec_lsize(), PETSC_DECIDE, v_rhs, &rhs);
57 if (sol_array == NULL) v_solution=(
double *)xmalloc(
sizeof(
double) * (this->
vec_lsize() + 1));
58 else v_solution=sol_array;
59 VecCreateMPIWithArray(PETSC_COMM_WORLD,1, this->
vec_lsize(), PETSC_DECIDE, v_solution, &solution);
72 finalize(MAT_FLUSH_ASSEMBLY);
78 ASSERT_PERMANENT(0).error(
"Can not set values. Matrix is not preallocated.\n");
90 finalize(MAT_FLUSH_ASSEMBLY);
96 ASSERT_PERMANENT(0).error(
"Can not set values. Matrix is not preallocated.\n");
102 void LinSys::finalize(MatAssemblyType assembly_type)
105 WarningOut() <<
"Finalizing linear system without setting values.\n";
106 preallocate_matrix();
108 MatAssemblyBegin(matrix, assembly_type);
109 VecAssemblyBegin(rhs);
110 if (assembly_type == MAT_FINAL_ASSEMBLY) status=
DONE;
111 MatAssemblyEnd(matrix, assembly_type);
115 void view(std::ostream output_stream,
int * output_mapping = NULL)
123 if (matrix)
chkerr(MatDestroy(&matrix));
125 VecDestroy(&solution);
128 if (own_solution) xfree(v_solution);
135 void LSView(LinSystem *ls)
140 const PetscInt *
cols;
141 const PetscScalar *vals;
145 MatGetInfo(ls->A,MAT_GLOBAL_SUM,&info);
147 if (ls->ds->myp() == 0) {
148 f=xfopen(
"matrix.dat",
"w");
149 xfprintf(f,
"%d %d 0.0\n",ls->size, ls->size);
154 for(n=0;n<ls->ds->np();n++) {
156 if (n==ls->ds->myp()) {
157 f=xfopen(
"matrix.dat",
"a");
158 for(i=ls->ds->begin(); i<ls->ds->end();i++) {
159 MatGetRow(ls->A,i,&ncols,&
cols,&vals);
161 xfprintf(f,
"%d %d %f\n",ls->old_4_new[i]+1,ls->old_4_new[
cols[j]]+1,vals[j]);
162 MatRestoreRow(ls->A,i,&ncols,&
cols,&vals);
168 if (ls->ds->myp() == 0) {
169 f=xfopen(
"rhs.m",
"w");
170 xfprintf(f,
"yyy = zeros(%d,2)\n yyy=[\n",ls->size);
174 for(n=0; n<ls->ds->np(); n++) {
176 if (n==ls->ds->myp()) {
177 f=xfopen(
"rhs.m",
"a");
178 for(i=ls->ds->begin(); i<ls->ds->end(); i++) {
179 xfprintf(f,
"%d %f\n",ls->old_4_new[i]+1,ls->vb[i-ls->ds->begin()]);
185 if (ls->ds->myp() == 0) {
186 f=xfopen(
"rhs.m",
"a");
187 xfprintf(f,
"];\n xxx=sortrows(yyy); Vec=xxx(:,2);\n");
195 void MyVecView(Vec v,
int *order,
const char* fname) {
209 xfprintf(f,
"yyy = zeros(%d,2)\n yyy=[\n",ds.size());
213 for(n=0; n<ds.np(); n++) {
216 VecGetArray(v,&
array);
218 for(i=ds.begin(); i<ds.end(); i++) {
219 xfprintf(f,
"%d %f\n",order[i]+1,
array[i-ds.begin()]);
221 VecRestoreArray(v,&
array);
228 xfprintf(f,
"];\n xxx=sortrows(yyy); Vec=xxx(:,2);\n");
242 void LSSetCSR( LinSystem *mtx )
246 const PetscInt *
cols;
247 const PetscScalar *vals;
250 MatGetInfo(mtx->A,MAT_LOCAL,&info);
251 nnz=(int)(info.nz_used)+1;
253 if (mtx->i == NULL) mtx->i=(
int *)xmalloc( (mtx->size+1)*
sizeof(int) );
254 if (mtx->j == NULL) mtx->j=(
int *)xmalloc( nnz*
sizeof(
int) );
255 if (mtx->a == NULL) mtx->a=(
double *)xmalloc( nnz*
sizeof(
double) );
256 DebugOut().fmt(
"mtx size: {} nnz: {}\n", mtx->size, nnz);
259 for( row=0; row< mtx->size; row++ ) {
261 MatGetRow(mtx->A,row,&nnz_loc,&
cols,&vals);
262 DebugOut().fmt(
"CSR row: {} nnz_loc {}\n", row, nnz_loc);
263 for(i=0; i<nnz_loc; i++) {
264 ASSERT_PERMANENT(pos<nnz)(row)(i).error(
"More nonzeroes then allocated!\n");
269 MatRestoreRow(mtx->A,row,&nnz_loc,&
cols,&vals);
274 void LSFreeCSR( LinSystem *mtx )
285 void LinSys_MPIAIJ::start_allocation()
287 if (status != NONE) {
291 VecCreateMPI(PETSC_COMM_WORLD, vec_ds.lsize(), PETSC_DECIDE, &(on_vec));
292 VecDuplicate(on_vec,&(off_vec));
297 void LinSys_MPIAIJ::preallocate_matrix()
299 ASSERT_PERMANENT(status == ALLOCATE).error(
"Linear system has to be in ALLOCATE status.");
301 PetscScalar *on_array, *off_array;
306 VecAssemblyBegin(on_vec);
307 VecAssemblyBegin(off_vec);
309 on_nz=(
int *)xmalloc( 2 * vec_ds.lsize() *
sizeof(int));
310 off_nz=on_nz + vec_ds.lsize();
312 VecAssemblyEnd(on_vec);
313 VecAssemblyEnd(off_vec);
315 VecGetArray(on_vec,&on_array);
316 VecGetArray(off_vec,&off_array);
318 for(i=0; i<vec_ds.lsize(); i++) {
319 on_nz[i]=min((
unsigned int)(on_array[i]+0.1),vec_ds.lsize());
320 off_nz[i]=(int)(off_array[i]+0.1);
323 VecRestoreArray(on_vec,&on_array);
324 VecRestoreArray(off_vec,&off_array);
326 VecDestroy(&off_vec);
329 MatCreateAIJ(PETSC_COMM_WORLD, vec_ds.lsize(), vec_ds.lsize(), PETSC_DETERMINE, PETSC_DETERMINE,
330 0, on_nz, 0, off_nz, &matrix);
332 if (symmetric) MatSetOption(matrix, MAT_SYMMETRIC, PETSC_TRUE);
333 MatSetOption(matrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
339 void LinSys_MPIAIJ::preallocate_values(
int nrow,
int *
rows,
int ncol,
int *
cols)
343 for (i=0; i<nrow; i++) {
345 for(j=0; j<ncol; j++) {
347 if (vec_ds.get_proc(row) == vec_ds.get_proc(col))
348 VecSetValue(on_vec,row,1.0,ADD_VALUES);
350 VecSetValue(off_vec,row,1.0,ADD_VALUES);
355 void LinSys_MPIAIJ::view_local_matrix()
358 MessageOut() <<
"Printing of local matrix is not supported yet for MPIAIJ matrix. \n";
363 LinSys_MPIAIJ:: ~LinSys_MPIAIJ()
365 if (status == ALLOCATE) {
367 VecDestroy(&off_vec);
373 LinSys_MATIS::LinSys_MATIS(std::shared_ptr<LocalToGlobalMap> global_row_4_sub_row,
double *sol_array)
374 :
LinSys(global_row_4_sub_row->get_distr()->lsize(), sol_array), lg_map(global_row_4_sub_row)
383 ASSERT_PERMANENT_LE(lg_map->get_distr()->size(), numeric_limits<PetscInt>::max()).error(
"Index range doesn't fit into signed int!");
384 err = ISLocalToGlobalMappingCreate(PETSC_COMM_WORLD,
386 (
const PetscInt*)(&(lg_map->get_map_vector()[0])),
387 PETSC_COPY_VALUES, &map_local_to_global);
393 loc_rows =
new int[loc_rows_size];
399 void LinSys_MATIS::start_allocation()
429 void LinSys_MATIS::preallocate_matrix()
431 ASSERT_PERMANENT(status == ALLOCATE).error(
"Linear system has to be in ALLOCATE status.");
438 MatSeqAIJSetPreallocation(local_matrix, 0, subdomain_nz);
440 if (symmetric) MatSetOption(matrix, MAT_SYMMETRIC, PETSC_TRUE);
441 MatSetOption(matrix, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
443 delete[] subdomain_nz;
446 void LinSys_MATIS::preallocate_values(
int nrow,
int *
rows,
int ncol,
int *
cols)
448 int i,row, n_loc_rows;
451 if (loc_rows_size < nrow) {
454 loc_rows=
new int[loc_rows_size];
459 err = ISGlobalToLocalMappingApply(map_local_to_global, IS_GTOLM_DROP, nrow,
rows, &n_loc_rows, loc_rows);
461 ASSERT_PERMANENT(nrow == n_loc_rows).error(
"Not all global indices translated to local indices.");
477 for (i=0; i<n_loc_rows; i++) {
479 subdomain_nz[row] = subdomain_nz[row] + ncol;
483 void LinSys_MATIS::view_local_matrix()
495 err = MatView(local_matrix,PETSC_VIEWER_STDOUT_SELF);
517 LinSys_MATIS:: ~LinSys_MATIS()
522 err = ISLocalToGlobalMappingDestroy(&map_local_to_global);
523 ASSERT_PERMANENT(err == 0).error(
"Error in ISLocalToGlobalMappingDestroy.");
528 delete[] subdomain_indices;
529 if (status == ALLOCATE) {
538 #ifdef HAVE_ATLAS_ONLY_LAPACK
545 void dgeqrf_(
int m,
int n,
double **A,
double *TAU,
double *work,
int lwork,
int info)