Commit f1382370 authored by kadir-hpi7's avatar kadir-hpi7

Trailing whitespaces are removed

parent 7f93d0b7
......@@ -7,14 +7,14 @@ csv zero = 0.0;
#define CS_MAX(a,b) (((a) > (b)) ? (a) : (b))
typedef struct csr_t {
csi m;
csi n;
csi nzmax;
csi nr;
csi* r;
csi* p;
csi* j;
csv* x;
csi m;
csi n;
csi nzmax;
csi nr;
csi* r;
csi* p;
csi* j;
csv* x;
} csr;
/* method declarations {{{*/
csr *csr_spfree(csr *A);
......@@ -24,14 +24,14 @@ csr *csr_spfree(csr *A);
/* free workspace and return a sparse matrix result */
csr *csr_done(csr *C, void *w, void *x, csi ok) {
return(ok ? C : csr_spfree(C)); /* return result if OK, else free it */
return (ok ? C : csr_spfree(C)); /* return result if OK, else free it */
}
/* wrapper for free */
void *cs_free(void *p) {
if (p)
free(p); /* free p if it is not already NULL */
return(NULL); /* return NULL to simplify the use of cs_free */
return (NULL); /* return NULL to simplify the use of cs_free */
}
/* wrapper for realloc */
......@@ -40,9 +40,9 @@ void *csr_realloc(void *p, csi n, size_t size, csi *ok) {
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("%d:reallocation failed, pnew is NULL\n", __LINE__);
printf("%d:reallocation failed, pnew is NULL\n", __LINE__);
}
return((*ok) ? pnew : p); /* return original p if failure */
return ((*ok) ? pnew : p); /* return original p if failure */
}
/* wrapper for realloc */
......@@ -51,31 +51,31 @@ void *cs_realloc(void *p, csi n, size_t size, csi *ok) {
pnew = realloc(p, CS_MAX(n, 1) * size); /* realloc the block */
*ok = (pnew != NULL); /* realloc fails if pnew is NULL */
if (pnew == NULL) {
printf("reallocation failed\n");
printf("reallocation failed\n");
}
return((*ok) ? pnew : p); /* return original p if failure */
return ((*ok) ? pnew : p); /* return original p if failure */
}
/* change the max # of entries sparse matrix */
csi csr_sprealloc(csr *A, csi nzmax) {
csi ok, oki = 0, okj = 1, okx = 1;
if (!A)
return(0);
return (0);
if (nzmax <= 0)
nzmax = A->p[A->m];
A->j = (int*)csr_realloc(A->j, nzmax, sizeof(csi), &oki);
nzmax = A->p[A->m];
A->j = (int*) csr_realloc(A->j, nzmax, sizeof (csi), &oki);
if (A->x)
A->x = (csv*)csr_realloc(A->x, nzmax, sizeof(csv), &okx);
A->x = (csv*) csr_realloc(A->x, nzmax, sizeof (csv), &okx);
ok = (oki && okj && okx);
if (ok)
A->nzmax = nzmax;
return(ok);
A->nzmax = nzmax;
return (ok);
}
/* free a sparse matrix */
csr *csr_spfree(csr *A) {
if (!A)
return(NULL); /* do nothing if A already NULL */
return (NULL); /* do nothing if A already NULL */
cs_free(A->p);
A->p = NULL;
cs_free(A->j);
......@@ -90,30 +90,30 @@ csr *csr_spfree(csr *A) {
/* allocate a sparse matrix (triplet form or compressed-ROW form) */
csr *csr_spalloc(csi m, csi n, csi nzmax, int values, int triplet, csv f) {
csr* A = (csr*)calloc(1, sizeof(csr)); /* allocate the cs struct */
csr* A = (csr*) calloc(1, sizeof (csr)); /* allocate the cs struct */
if (!A) {
perror("sparse allocation failed");
return(NULL); /* out of memory */
perror("sparse allocation failed");
return (NULL); /* out of memory */
}
A->m = m; /* define dimensions and nzmax */
A->n = n;
A->nzmax = nzmax = CS_MAX(nzmax, 0);
A->nr = 0; // number of nonzero rows
A->p = (csi*)calloc(m + 2, sizeof(csi));
A->j = (csi*)calloc(CS_MAX(nzmax,1), sizeof(csi));
A->x = (csv*)calloc(CS_MAX(nzmax,1), sizeof(csv));
return((!A->p || !A->j || !A->x) ? csr_spfree(A) : A);
A->nr = 0; // number of nonzero rows
A->p = (csi*) calloc(m + 2, sizeof (csi));
A->j = (csi*) calloc(CS_MAX(nzmax, 1), sizeof (csi));
A->x = (csv*) calloc(CS_MAX(nzmax, 1), sizeof (csv));
return ((!A->p || !A->j || !A->x) ? csr_spfree(A) : A);
}/*}}}*/
/** Multiply two sparse matrices which are stored in CSR format. MKL is used */
csr *csr_multiply(csi Am, csi An, csi Anzmax, const csi* Ap, const csi* Aj, const csv* Ax, csi Bm, csi Bn, csi Bnzmax, const csi* Bp, const csi* Bj, const csv* Bx, long* nummult, csi* xb, csv* x) { /*{{{*/
csv tf = 0;
csi p, jp, j, kp, k, i, nz = 0, anz, *Cp, *Cj, m, n,
bnz, values = 1;
bnz, values = 1;
csv *Cx;
csr *C;
if (An != Bm)
return(NULL);
return (NULL);
if (Anzmax == 0 || Bnzmax == 0) {
C = csr_spalloc(Am, Bn, 0, values, 0, tf);
return C;
......@@ -122,8 +122,8 @@ csr *csr_multiply(csi Am, csi An, csi Anzmax, const csi* Ap, const csi* Aj, cons
anz = Ap[Am];
n = Bn;
bnz = Bp[Bm];
for(i = 0; i < n; i++) xb[i] = 0;
for(i = 0; i < n; i++)
for (i = 0; i < n; i++) xb[i] = 0;
for (i = 0; i < n; i++)
xb[i] = 0;
values = (Ax != NULL) && (Bx != NULL);
csi tnz = (anz + bnz) * 2;
......@@ -132,8 +132,8 @@ csr *csr_multiply(csi Am, csi An, csi Anzmax, const csi* Ap, const csi* Aj, cons
return (csr_done(C, xb, x, 0));
Cp = C->p;
for (i = 0; i < m; i++) {
if ( ( (nz + n) > C->nzmax ) ) {
if(!csr_sprealloc(C, (2 * (C->nzmax) + n) ) ) {
if (((nz + n) > C->nzmax)) {
if (!csr_sprealloc(C, (2 * (C->nzmax) + n))) {
return (csr_done(C, xb, x, 0)); // out of memory
} else {
}
......
This diff is collapsed.
......@@ -8,28 +8,28 @@ MIC_TARGET int option_print_matrices = 0;
/** Converts COO matrix to CSR matrix */
void coo_to_csr(int m, int nnz, int* I, int* J, double* val, MKL_INT* AI, MKL_INT* AJ, double* Aval) { /*{{{*/
MKL_INT info = 0;
MKL_INT job[8];
job[1]=1; // one based indexing in csr
job[2]=0; // zero based indexing in coo
job[3]=2; // I don't know
job[4]=nnz; // nnz
job[0]=1; // coo to csr
job[5]=0; // Acsr and AJR allocated by user
mkl_dcsrcoo (job,&m, Aval, AJ, AI, &nnz, val, I, J, &info);
MKL_INT job[8];
job[1] = 1; // one based indexing in csr
job[2] = 0; // zero based indexing in coo
job[3] = 2; // I don't know
job[4] = nnz; // nnz
job[0] = 1; // coo to csr
job[5] = 0; // Acsr and AJR allocated by user
mkl_dcsrcoo(job, &m, Aval, AJ, AI, &nnz, val, I, J, &info);
} /* ENDOF coo_to_csr }}}*/
/** Prints matrix in CSR format */
void MIC_TARGET printmm_one(int m, double* Aval, int* AJ, int* AI){ //{{{
void MIC_TARGET printmm_one(int m, double* Aval, int* AJ, int* AI) { //{{{
if (option_print_matrices == OPTION_NOPRINT_MATRICES)
return;
int i;
for(i = 0; i < m; i++) {
printf("%d: ", i+1);
for (i = 0; i < m; i++) {
printf("%d: ", i + 1);
int j;
for(j = AI[i]-1; j < AI[i+1]-1; j++) {
for (j = AI[i] - 1; j < AI[i + 1] - 1; j++) {
printf("%d:%g ", AJ[j], Aval[j]);
}
printf("\n");
......@@ -38,20 +38,20 @@ void MIC_TARGET printmm_one(int m, double* Aval, int* AJ, int* AI){ //{{{
}//}}}
/** Writes matrix in CSR format in to a file using Matrix Market format */
void MIC_TARGET printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI){//{{{
void MIC_TARGET printfilemm_one(char* file, int m, int n, double* Aval, int* AJ, int* AI) {//{{{
FILE* f = fopen(file, "w");
if(f == NULL){
if (f == NULL) {
printf("%s %s %d: %s cannot be opened to write matrix\n", __FILE__, __PRETTY_FUNCTION__, __LINE__, file);
exit(1);
exit(1);
}
int i;
fprintf(f, "%%%%MatrixMarket matrix coordinate real general\n");
fprintf(f, "%d %d %d\n", m, n, AI[m]-1);
for(i = 0; i < m; i++) {
fprintf(f, "%d %d %d\n", m, n, AI[m] - 1);
for (i = 0; i < m; i++) {
int j;
for(j = AI[i]-1; j < AI[i+1]-1; j++) {
fprintf(f, "%d %d %g\n", i+1, AJ[j], Aval[j]);
for (j = AI[i] - 1; j < AI[i + 1] - 1; j++) {
fprintf(f, "%d %d %g\n", i + 1, AJ[j], Aval[j]);
}
}
fclose(f);
......
......@@ -7,9 +7,9 @@ All rights reserved.
Redistribution and use in source and binary forms, with or without modification,
are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright notice, this
* Redistributions of source code must retain the above copyright notice, this
list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright notice, this
* Redistributions in binary form must reproduce the above copyright notice, this
list of conditions and the following disclaimer in the documentation and/or
other materials provided with the distribution.
......@@ -23,7 +23,7 @@ are permitted provided that the following conditions are met:
ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
*/
*/
static char help[] = "Solves a sequence of linear systems with different right-hand-side vectors.\n\
Input parameters include:\n\
......@@ -46,176 +46,232 @@ T*/
petscmat.h - matrices
petscis.h - index sets petscksp.h - Krylov subspace methods
petscviewer.h - viewers petscpc.h - preconditioners
*/
*/
#include <petscksp.h>
#undef __FUNCT__
#define __FUNCT__ "main"
int main(int argc,char **args)
{
Vec x,b,u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* linear solver context */
PetscReal norm; /* norm of solution error */
PetscErrorCode ierr;
PetscInt ntimes,i,j,k,Ii,J,Istart,Iend;
PetscInt m = 8,n = 7,its;
PetscBool flg = PETSC_FALSE;
PetscScalar v,one = 1.0,neg_one = -1.0,rhs;
PetscInitialize(&argc,&args,(char*)0,help);
ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrix for use in solving a series of
linear systems of the form, A x_i = b_i, for i=1,2,...
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create parallel matrix, specifying only its global dimensions.
When using MatCreate(), the matrix format can be specified at
runtime. Also, the parallel partitioning of the matrix is
determined by PETSc at runtime.
*/
ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr);
ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr);
ierr = MatSetFromOptions(A);CHKERRQ(ierr);
ierr = MatSetUp(A);CHKERRQ(ierr);
/*
Currently, all PETSc parallel matrix formats are partitioned by
contiguous chunks of rows across the processors. Determine which
rows of the matrix are locally owned.
*/
ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr);
/*
Set matrix elements for the 2-D, five-point stencil in parallel.
- Each processor needs to insert only elements that it owns
locally (but any non-local elements will be sent to the
appropriate processor during matrix assembly).
- Always specify global rows and columns of matrix entries.
*/
for (Ii=Istart; Ii<Iend; Ii++) {
v = -1.0; i = Ii/n; j = Ii - i*n;
if (i>0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (i<m-1) {J = Ii + n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j>0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
if (j<n-1) {J = Ii + 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);}
v = 4.0; ierr = MatSetValues(A,1,&Ii,1,&Ii,&v,INSERT_VALUES);CHKERRQ(ierr);
}
/*
Assemble matrix, using the 2-step process:
MatAssemblyBegin(), MatAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
/*
Create parallel vectors.
- When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
we specify only the vector's global
dimension; the parallel partitioning is determined at runtime.
- When solving a linear system, the vectors and matrices MUST
be partitioned accordingly. PETSc automatically generates
appropriately partitioned matrices and vectors when MatCreate()
and VecCreate() are used with the same communicator.
- Note: We form 1 vector from scratch and then duplicate as needed.
*/
ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr);
ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr);
ierr = VecSetFromOptions(u);CHKERRQ(ierr);
ierr = VecDuplicate(u,&b);CHKERRQ(ierr);
ierr = VecDuplicate(b,&x);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);
/*
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);
/*
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve several linear systems of the form A x_i = b_i
I.e., we retain the same matrix (A) for all systems, but
change the right-hand-side vector (b_i) at each step.
In this case, we simply call KSPSolve() multiple times. The
preconditioner setup operations (e.g., factorization for ILU)
be done during the first call to KSPSolve() only; such operations
will NOT be repeated for successive solves.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ntimes = 2;
ierr = PetscOptionsGetInt(NULL,"-ntimes",&ntimes,NULL);CHKERRQ(ierr);
for (k=1; k<ntimes+1; k++) {
int main(int argc, char **args) {
Vec x, b, u; /* approx solution, RHS, exact solution */
Mat A; /* linear system matrix */
KSP ksp; /* linear solver context */
PetscReal norm; /* norm of solution error */
PetscErrorCode ierr;
PetscInt ntimes, i, j, k, Ii, J, Istart, Iend;
PetscInt m = 8, n = 7, its;
PetscBool flg = PETSC_FALSE;
PetscScalar v, one = 1.0, neg_one = -1.0, rhs;
PetscInitialize(&argc, &args, (char*) 0, help);
ierr = PetscOptionsGetInt(NULL, "-m", &m, NULL);
CHKERRQ(ierr);
ierr = PetscOptionsGetInt(NULL, "-n", &n, NULL);
CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Compute the matrix for use in solving a series of
linear systems of the form, A x_i = b_i, for i=1,2,...
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create parallel matrix, specifying only its global dimensions.
When using MatCreate(), the matrix format can be specified at
runtime. Also, the parallel partitioning of the matrix is
determined by PETSc at runtime.
*/
ierr = MatCreate(PETSC_COMM_WORLD, &A);
CHKERRQ(ierr);
ierr = MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, m*n, m * n);
CHKERRQ(ierr);
ierr = MatSetFromOptions(A);
CHKERRQ(ierr);
ierr = MatSetUp(A);
CHKERRQ(ierr);
/*
Currently, all PETSc parallel matrix formats are partitioned by
contiguous chunks of rows across the processors. Determine which
rows of the matrix are locally owned.
*/
ierr = MatGetOwnershipRange(A, &Istart, &Iend);
CHKERRQ(ierr);
/*
Set matrix elements for the 2-D, five-point stencil in parallel.
- Each processor needs to insert only elements that it owns
locally (but any non-local elements will be sent to the
appropriate processor during matrix assembly).
- Always specify global rows and columns of matrix entries.
*/
for (Ii = Istart; Ii < Iend; Ii++) {
v = -1.0;
i = Ii / n;
j = Ii - i*n;
if (i > 0) {
J = Ii - n;
ierr = MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
CHKERRQ(ierr);
}
if (i < m - 1) {
J = Ii + n;
ierr = MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
CHKERRQ(ierr);
}
if (j > 0) {
J = Ii - 1;
ierr = MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
CHKERRQ(ierr);
}
if (j < n - 1) {
J = Ii + 1;
ierr = MatSetValues(A, 1, &Ii, 1, &J, &v, INSERT_VALUES);
CHKERRQ(ierr);
}
v = 4.0;
ierr = MatSetValues(A, 1, &Ii, 1, &Ii, &v, INSERT_VALUES);
CHKERRQ(ierr);
}
/*
Assemble matrix, using the 2-step process:
MatAssemblyBegin(), MatAssemblyEnd()
Computations can be done while messages are in transition
by placing code between these two statements.
*/
ierr = MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
ierr = MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
CHKERRQ(ierr);
/*
Create parallel vectors.
- When using VecCreate(), VecSetSizes() and VecSetFromOptions(),
we specify only the vector's global
dimension; the parallel partitioning is determined at runtime.
- When solving a linear system, the vectors and matrices MUST
be partitioned accordingly. PETSc automatically generates
appropriately partitioned matrices and vectors when MatCreate()
and VecCreate() are used with the same communicator.
- Note: We form 1 vector from scratch and then duplicate as needed.
*/
ierr = VecCreate(PETSC_COMM_WORLD, &u);
CHKERRQ(ierr);
ierr = VecSetSizes(u, PETSC_DECIDE, m * n);
CHKERRQ(ierr);
ierr = VecSetFromOptions(u);
CHKERRQ(ierr);
ierr = VecDuplicate(u, &b);
CHKERRQ(ierr);
ierr = VecDuplicate(b, &x);
CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Create the linear solver and set various options
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Create linear solver context
*/
ierr = KSPCreate(PETSC_COMM_WORLD, &ksp);
CHKERRQ(ierr);
/*
Set exact solution; then compute right-hand-side vector. We use
an exact solution of a vector with all elements equal to 1.0*k.
*/
rhs = one * (PetscReal)k;
ierr = VecSet(u,rhs);CHKERRQ(ierr);
ierr = MatMult(A,u,b);CHKERRQ(ierr);
Set operators. Here the matrix that defines the linear system
also serves as the preconditioning matrix.
*/
ierr = KSPSetOperators(ksp, A, A);
CHKERRQ(ierr);
/*
View the exact solution vector if desired
*/
ierr = PetscOptionsGetBool(NULL,"-view_exact_sol",&flg,NULL);CHKERRQ(ierr);
if (flg) {ierr = VecView(u,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);}
Set runtime options, e.g.,
-ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
These options will override those specified above as long as
KSPSetFromOptions() is called _after_ any other customization
routines.
*/
ierr = KSPSetFromOptions(ksp);
CHKERRQ(ierr);
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Solve several linear systems of the form A x_i = b_i
I.e., we retain the same matrix (A) for all systems, but
change the right-hand-side vector (b_i) at each step.
In this case, we simply call KSPSolve() multiple times. The
preconditioner setup operations (e.g., factorization for ILU)
be done during the first call to KSPSolve() only; such operations
will NOT be repeated for successive solves.
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr);
ntimes = 2;
ierr = PetscOptionsGetInt(NULL, "-ntimes", &ntimes, NULL);
CHKERRQ(ierr);
for (k = 1; k < ntimes + 1; k++) {
/*
Set exact solution; then compute right-hand-side vector. We use
an exact solution of a vector with all elements equal to 1.0*k.
*/
rhs = one * (PetscReal) k;
ierr = VecSet(u, rhs);
CHKERRQ(ierr);
ierr = MatMult(A, u, b);
CHKERRQ(ierr);
/*
View the exact solution vector if desired
*/
ierr = PetscOptionsGetBool(NULL, "-view_exact_sol", &flg, NULL);
CHKERRQ(ierr);
if (flg) {
ierr = VecView(u, PETSC_VIEWER_STDOUT_WORLD);
CHKERRQ(ierr);
}
ierr = KSPSolve(ksp, b, x);
CHKERRQ(ierr);
/*
Check the error
*/
ierr = VecAXPY(x, neg_one, u);
CHKERRQ(ierr);
ierr = VecNorm(x, NORM_2, &norm);
CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp, &its);
CHKERRQ(ierr);
/*
Print convergence information. PetscPrintf() produces a single
print statement from all processes that share a communicator.
*/
ierr = PetscPrintf(PETSC_COMM_WORLD, "Norm of error %g System %D: iterations %D\n", (double) norm, k, its);
CHKERRQ(ierr);
}
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Clean up
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
/*
Check the error
*/
ierr = VecAXPY(x,neg_one,u);CHKERRQ(ierr);
ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr);
ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr);
Free work space. All PETSc objects should be destroyed when they
are no longer needed.
*/
ierr = KSPDestroy(&ksp);
CHKERRQ(ierr);
ierr = VecDestroy(&u);
CHKERRQ(ierr);