Actual source code: iccbs.c

      1: #define PETSCMAT_DLL

      3: /*
      4:    Defines a Cholesky factorization preconditioner with BlockSolve95 interface.

      6:    Note that BlockSolve95 works with a scaled and permuted preconditioning matrix.
      7:    If the linear system matrix and preconditioning matrix are the same, we then
      8:    work directly with the permuted and scaled linear system:
      9:       - original system:  Ax = b
     10:       - permuted and scaled system:   Cz = f, where
     11:              C = P D^{-1/2} A D^{-1/2}
     12:              z = P D^{1/2} x
     13:              f = P D^{-1/2} b
     14:              D = diagonal of A
     15:              P = permutation matrix determined by coloring
     16:    In this case, we use pre-solve and post-solve phases to handle scaling and
     17:    permutation, and by default the scaled residual norm is monitored for the
     18:    ILU/ICC preconditioners.  Use the option
     19:      -ksp_truemonitor
     20:    to print both the scaled and unscaled residual norms.
     21: */

     23:  #include petsc.h

     25:  #include src/mat/impls/rowbs/mpi/mpirowbs.h

     29: PetscErrorCode MatScaleSystem_MPIRowbs(Mat mat,Vec x,Vec rhs)
     30: {
     31:   Mat_MPIRowbs *bsif  = (Mat_MPIRowbs*)mat->data;
     32:   Vec          v = bsif->xwork;
     33:   PetscScalar  *xa,*rhsa,*va;

     37:   /* Permute and scale RHS and solution vectors */
     38:   if (x) {
     39:     VecGetArray(x,&xa);
     40:     VecGetArray(v,&va);
     41:     BSperm_dvec(xa,va,bsif->pA->perm);CHKERRBS(0);
     42:     VecRestoreArray(x,&xa);
     43:     VecRestoreArray(v,&va);
     44:     VecPointwiseDivide(x,v,bsif->diag);
     45:   }

     47:   if (rhs) {
     48:     VecGetArray(rhs,&rhsa);
     49:     VecGetArray(v,&va);
     50:     BSperm_dvec(rhsa,va,bsif->pA->perm);CHKERRBS(0);
     51:     VecRestoreArray(rhs,&rhsa);
     52:     VecRestoreArray(v,&va);
     53:     VecPointwiseMult(rhs,v,bsif->diag);
     54:   }
     55:   return(0);
     56: }

     60: PetscErrorCode MatUnScaleSystem_MPIRowbs(Mat mat,Vec x,Vec rhs)
     61: {
     62:   Mat_MPIRowbs *bsif  = (Mat_MPIRowbs*)mat->data;
     63:   Vec          v = bsif->xwork;
     64:   PetscScalar  *xa,*va,*rhsa;

     68:   /* Unpermute and unscale the solution and RHS vectors */
     69:   if (x) {
     70:     VecPointwiseMult(v,x,bsif->diag);
     71:     VecGetArray(v,&va);
     72:     VecGetArray(x,&xa);
     73:     BSiperm_dvec(va,xa,bsif->pA->perm);CHKERRBS(0);
     74:     VecRestoreArray(x,&xa);
     75:     VecRestoreArray(v,&va);
     76:   }
     77:   if (rhs) {
     78:     VecPointwiseDivide(v,rhs,bsif->diag);
     79:     VecGetArray(rhs,&rhsa);
     80:     VecGetArray(v,&va);
     81:     BSiperm_dvec(va,rhsa,bsif->pA->perm);CHKERRBS(0);
     82:     VecRestoreArray(rhs,&rhsa);
     83:     VecRestoreArray(v,&va);
     84:   }
     85:   return(0);
     86: }

     90: PetscErrorCode MatUseScaledForm_MPIRowbs(Mat mat,PetscTruth scale)
     91: {
     92:   Mat_MPIRowbs *bsif  = (Mat_MPIRowbs*)mat->data;

     95:   bsif->vecs_permscale = scale;
     96:   return(0);
     97: }