<html><head><meta http-equiv="content-type" content="text/html; charset=us-ascii"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br></div>It is pasted below<br><div><br><blockquote type="cite"><div>On Jun 25, 2023, at 1:40 PM, Edoardo alinovi <edoardo.alinovi@gmail.com> wrote:</div><br class="Apple-interchange-newline"><div><div dir="ltr">Hi Barry, <div><br></div><div>thanks for pointing me out to that discussion! Unfortunately I am getting issues with this link:
<a href="http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190302/b0c1ad29/attachment.mht" target="_blank" style="white-space-collapse: preserve;">http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20190302/b0c1ad29/attachment.mht</a> , any chance it is a dead one?<div><br></div><div>Cheers</div></div></div>
</div></blockquote></div><br><div><span style="font-size: 13px;">--====-=-=</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Type: text/plain</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Disposition: inline</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> Ok, I have implemented the algorithm using PETSc PCCOMPOSITE and PCGALERKIN and get identical iterations as the code you sent. PCFIELDSPLIT is not intended for this type of solver composition. Here is the algorithm written in "two-step" form</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> x_1/2 = P KSPSolve( R A P, using BoomerAMG) R b</span><br style="font-size: 13px;"><span style="font-size: 13px;"> x = x_1/2 + PCApply( A, using Hypre PILUT preconditioner) ( b - A x_1/2)</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">PCCOMPOSITE with a type of multiplicative handles the two steps and PCGALERKIN handles the P KSPSolve(R A P) R business.</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">You will need to use the master version of PETSc because I had to add a feature to PCGALERKIN to allow the solver to be easily used for and A that changes values for later solvers.</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">Here is the output from -ksp_view </span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">KSP Object: 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: fgmres</span><br style="font-size: 13px;"><span style="font-size: 13px;"> GMRES: restart=100, using Classical (unmodified) Gram-Schmidt Orthogonalization with no iterative refinement</span><br style="font-size: 13px;"><span style="font-size: 13px;"> GMRES: happy breakdown tolerance 1e-30</span><br style="font-size: 13px;"><span style="font-size: 13px;"> maximum iterations=10000, initial guess is zero</span><br style="font-size: 13px;"><span style="font-size: 13px;"> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> right preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;"> using UNPRECONDITIONED norm type for convergence test</span><br style="font-size: 13px;"><span style="font-size: 13px;">PC Object: 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: composite</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Composite PC type - MULTIPLICATIVE</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PCs on composite preconditioner follow</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ---------------------------------</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC Object: (sub_0_) 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: galerkin</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Galerkin PC</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSP on Galerkin follow</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ---------------------------------</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSP Object: (sub_0_galerkin_) 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: richardson</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Richardson: damping factor=1.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> maximum iterations=10000, initial guess is zero</span><br style="font-size: 13px;"><span style="font-size: 13px;"> tolerances: relative=1e-05, absolute=1e-50, divergence=10000.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> left preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;"> using PRECONDITIONED norm type for convergence test</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC Object: (sub_0_galerkin_) 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: hypre</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Cycle type V</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Maximum number of levels 25</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Maximum number of iterations PER hypre call 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Convergence tolerance PER hypre call 0.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Threshold for strong coupling 0.25</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Interpolation truncation factor 0.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Interpolation: max elements per row 0</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Number of levels of aggressive coarsening 0</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Number of paths for aggressive coarsening 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Maximum row sums 0.9</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Sweeps down 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Sweeps up 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Sweeps on coarse 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Relax down symmetric-SOR/Jacobi</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Relax up symmetric-SOR/Jacobi</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Relax on coarse Gaussian-elimination</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Relax weight (all) 1.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Outer relax weight (all) 1.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Using CF-relaxation</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Not using more complex smoothers.</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Measure type local</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Coarsen type Falgout</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE BoomerAMG: Interpolation type classical</span><br style="font-size: 13px;"><span style="font-size: 13px;"> linear system matrix = precond matrix:</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat Object: 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: seqaij</span><br style="font-size: 13px;"><span style="font-size: 13px;"> rows=50, cols=50</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total: nonzeros=244, allocated nonzeros=244</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total number of mallocs used during MatSetValues calls =0</span><br style="font-size: 13px;"><span style="font-size: 13px;"> not using I-node routines</span><br style="font-size: 13px;"><span style="font-size: 13px;"> linear system matrix = precond matrix:</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat Object: 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: seqaij</span><br style="font-size: 13px;"><span style="font-size: 13px;"> rows=100, cols=100, bs=2</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total: nonzeros=976, allocated nonzeros=100000</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total number of mallocs used during MatSetValues calls =0</span><br style="font-size: 13px;"><span style="font-size: 13px;"> using I-node routines: found 50 nodes, limit used is 5</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC Object: (sub_1_) 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: hypre</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE Pilut preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE Pilut: maximum number of iterations 1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE Pilut: drop tolerance 0.1</span><br style="font-size: 13px;"><span style="font-size: 13px;"> HYPRE Pilut: default factor row size </span><br style="font-size: 13px;"><span style="font-size: 13px;"> linear system matrix = precond matrix:</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat Object: 1 MPI processes</span><br style="font-size: 13px;"><span style="font-size: 13px;"> type: seqaij</span><br style="font-size: 13px;"><span style="font-size: 13px;"> rows=100, cols=100, bs=2</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total: nonzeros=976, allocated nonzeros=100000</span><br style="font-size: 13px;"><span style="font-size: 13px;"> total number of mallocs used during MatSetValues calls =0</span><br style="font-size: 13px;"><span style="font-size: 13px;"> using I-node routines: found 50 nodes, limit used is 5</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">Note that one can change the solvers on the two "stages" and all their options such as tolerances etc using the options database a proper prefixes which you can find in the output above. For example to use PETSc's ILU instead of hypre's just run with -sub_1_pc_type ilu or to use a direct solver instead of boomer -sub_0_galerkin_pc_type lu </span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">I've attached a new version of testmain2.c that runs your solver and also my version.</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">Given the "unique" nature of your R = [ I I ... I] and P = [0 I 0 0 ... 0] I am not sure that it makes sense to include this preconditioner directly in PETSc as a new PC type; so if you are serious about using it you can take what I am sending back and modify it for your needs. As a numerical analyst who works on linear solvers I am also not convinced that this is likely to be a particular good preconditioner</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> Let me know if you have any questions,</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> Barry</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">--====-=-=</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Type: text/plain; name=testmain2.c</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Disposition: attachment; filename=testmain2.c</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">// make && mpirun -np 2 ./testmain2 -ksp_error_if_not_converged</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">#include "MCPR.h"</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">/*</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Computes the submatrix associated with the Galerkin subproblem Ap = R A P</span><br style="font-size: 13px;"><span style="font-size: 13px;">*/</span><br style="font-size: 13px;"><span style="font-size: 13px;">PetscErrorCode ComputeSubmatrix(PC pc,Mat A, Mat Ap, Mat *cAp,void *ctx)</span><br style="font-size: 13px;"><span style="font-size: 13px;">{</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscErrorCode ierr;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscInt b,Am,An,start,end,first = 0, offset = 1;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> IS is,js;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat Aij;</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionBegin;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetLocalSize(A,&Am,&An);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetBlockSize(A,&b);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetOwnershipRange(A, &start, &end);CHKERRQ(ierr);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISCreateStride(PetscObjectComm((PetscObject)A),Am/b,start+offset,b,&js);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> if (!Ap) {</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISCreateStride(PetscObjectComm((PetscObject)A),An/b,start+0,b,&is);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetSubMatrix(A,is,js,MAT_INITIAL_MATRIX,&Ap);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISDestroy(&is);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> *cAp = Ap;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> first = 1;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> } else {</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatZeroEntries(Ap);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> }</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> for(PetscInt k=first;k<b;++k) {</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISCreateStride(PetscObjectComm((PetscObject)A),An/b,start+k,b,&is);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetSubMatrix(A,is,js,MAT_INITIAL_MATRIX,&Aij);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatAXPY(Ap,1.0,Aij,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatDestroy(&Aij);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISDestroy(&is);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> }</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = ISDestroy(&js);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionReturn(0);</span><br style="font-size: 13px;"><span style="font-size: 13px;">}</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">/*</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Apply the restriction operator for the Galkerin problem</span><br style="font-size: 13px;"><span style="font-size: 13px;">*/</span><br style="font-size: 13px;"><span style="font-size: 13px;">PetscErrorCode ApplyR(Mat A, Vec x,Vec y)</span><br style="font-size: 13px;"><span style="font-size: 13px;">{</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscErrorCode ierr;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscInt b;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionBegin;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = VecGetBlockSize(x,&b);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = VecStrideGather(x,0,y,INSERT_VALUES);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> for (PetscInt k=1;k<b;++k) {ierr = VecStrideGather(x,k,y,ADD_VALUES);CHKERRQ(ierr);}</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionReturn(0);</span><br style="font-size: 13px;"><span style="font-size: 13px;">}</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">/*</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Apply the interpolation operator for the Galerkin problem</span><br style="font-size: 13px;"><span style="font-size: 13px;">*/</span><br style="font-size: 13px;"><span style="font-size: 13px;">PetscErrorCode ApplyP(Mat A, Vec x,Vec y)</span><br style="font-size: 13px;"><span style="font-size: 13px;">{</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscErrorCode ierr;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscInt offset = 1;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionBegin;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = VecStrideScatter(x,offset,y,INSERT_VALUES);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscFunctionReturn(0);</span><br style="font-size: 13px;"><span style="font-size: 13px;">}</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">int main( int argc, char *argv[] )</span><br style="font-size: 13px;"><span style="font-size: 13px;">{</span><br style="font-size: 13px;"><span style="font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscInitialize(&argc,&argv,PETSC_NULL,PETSC_NULL);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">int rank, size;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">MPI_Comm_rank (MPI_COMM_WORLD, &rank);</span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">/* get current process id */</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">MPI_Comm_size (MPI_COMM_WORLD, &size);</span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">/* get number of processes */</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">MPI_Comm C = PETSC_COMM_WORLD;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscRandom rnd;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscRandomCreate(C,&rnd);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscRandomSetInterval(rnd,0.0,1.0);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscRandomSetFromOptions(rnd);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">int M = 100;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">int N = size*M;</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Mat A = getSampleMatrix(M);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Mat T1 = getT1(A,1);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Mat T2 = getT2(A,0.1);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Mat MCPR = getMCPR(A,T1,T2);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Vec u,v,w,z;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecCreate(C,&u);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecSetBlockSize(u,2);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecSetSizes(u,M,N);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecSetFromOptions(u);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecDuplicate(u,&v);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecDuplicate(u,&w);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecDuplicate(u,&z);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">VecSetRandom(u,rnd);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">Mat mats[] = {T2,MCPR};</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">const char *names[] = {"T2","MCPR"};</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">for(int k=1;k<2;++k) {</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSP solver;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPCreate(C,&solver);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPSetOperators(solver,A,A);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPSetType(solver,KSPFGMRES);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPGMRESSetRestart(solver,N);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PC pc;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPGetPC(solver,&pc);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">putMatrixInPC(pc,mats[k]);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSPSetFromOptions(solver);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSPSolve(solver,u,v);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPConvergedReason reason;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">int its;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPGetConvergedReason(solver,&reason);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPGetIterationNumber(solver,&its);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">PetscPrintf(PETSC_COMM_WORLD,"testmain2: %s converged reason %d; iterations %d.\n",names[k],reason,its);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;"> KSPView(solver,PETSC_VIEWER_STDOUT_WORLD);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">KSPDestroy(&solver);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">}</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> /*</span><br style="font-size: 13px;"><span style="font-size: 13px;"> The preconditioner is</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> x_1/2 = P KSPSolve( R A P, using BoomerAMG) R b</span><br style="font-size: 13px;"><span style="font-size: 13px;"> x = x_1/2 + PCApply( A, using Hypre PILUT preconditioner) ( b - A x_1/2)</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> where the first line is implemented using PCGALERKIN with BoomerAMG on the subproblem so can be written as</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> x_1/2 = PCApply(A, using PCGALERKIN with KSPSolve( R A P, using BoomerAMG) as the inner solver</span><br style="font-size: 13px;"><span style="font-size: 13px;"> x = x_1/2 + PCApply( A, using Hypre PILUT preconditioner) ( b - A x_1/2)</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> Which is implemented using the PETSc PCCOMPOSITE preconditioner of type multiplicative so can be written as</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> x = PCApply(A, using PCCOMPOSITE using (PCGALERKIN with KSPSolve( R A P, using BoomerAMG) as the inner solver) as the first solver and PCApply( A, using Hypre PILUT preconditioner) as the second solver)</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> {</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscErrorCode ierr;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSP ksp;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPCreate(PETSC_COMM_WORLD,&ksp);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPSetType(ksp,KSPFGMRES);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPGMRESSetRestart(ksp,100);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPSetOperators(ksp,A,A);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC pc;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPGetPC(ksp,&pc);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCSetType(pc,PCCOMPOSITE);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCCompositeSetType(pc,PC_COMPOSITE_MULTIPLICATIVE);CHKERRQ(ierr);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> /* Create first sub problem solver Hypre boomerAMG on Ap */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC t1;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCCompositeAddPC(pc,PCGALERKIN);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCCompositeGetPC(pc,0,&t1);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> KSP Ap_ksp;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCGalerkinGetKSP(t1,&Ap_ksp);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPSetType(Ap_ksp,KSPRICHARDSON);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC Ap_pc;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPGetPC(Ap_ksp,&Ap_pc);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCSetType(Ap_pc,PCHYPRE);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> /* this tells the PC how to compute the reduced matrix */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCGalerkinSetComputeSubmatrix(t1,ComputeSubmatrix,NULL);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PetscInt b,Am,An;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = MatGetLocalSize(A,&Am,&An);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatGetBlockSize(A,&b);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">int start,end;</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = MatGetOwnershipRange(A, &start, &end);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> /* create the R operator */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat R;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatCreateShell(PetscObjectComm((PetscObject)A),Am/b,An,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&R);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatShellSetOperation(R,MATOP_MULT,(void (*)(void))ApplyR);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCGalerkinSetRestriction(t1,R);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> /* create the P operator */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> Mat P;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatCreateShell(PetscObjectComm((PetscObject)A),Am,An/b,PETSC_DETERMINE,PETSC_DETERMINE,NULL,&P);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = MatShellSetOperation(P,MATOP_MULT,(void (*)(void))ApplyP);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCGalerkinSetInterpolation(t1,P);CHKERRQ(ierr);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> /* Create the second subproblem solver Block ILU */</span><br style="font-size: 13px;"><span style="font-size: 13px;"> PC t2;</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCCompositeAddPC(pc,PCHYPRE);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = PCCompositeGetPC(pc,1,&t2);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = PCSetType(t2,PCHYPRE);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = PetscOptionsSetValue(NULL,"-sub_1_pc_hypre_pilut_maxiter","1");</span><br style="font-size: 13px;"><span style="font-size: 13px;"> char s[100];</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">sprintf(s,"%e",.1);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = PetscOptionsSetValue(NULL,"-sub_1_pc_hypre_pilut_tol",s);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">ierr = PCHYPRESetType(t2,"pilut");CHKERRQ(ierr);</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = VecZeroEntries(v);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPSolve(ksp,u,v);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> ierr = KSPView(ksp,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);</span><br style="font-size: 13px;"><span style="font-size: 13px;"> }</span><br style="font-size: 13px;"><span class="Apple-tab-span" style="white-space: pre; font-size: 13px;"> </span><span style="font-size: 13px;">return 0;</span><br style="font-size: 13px;"><span style="font-size: 13px;">}</span><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">--====-=-=</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Type: text/plain; charset=utf-8</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Disposition: inline</span><br style="font-size: 13px;"><span style="font-size: 13px;">Content-Transfer-Encoding: quoted-printable</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">> On Jan 23, 2017, at 8:25 AM, S=C3=A9bastien Loisel <</span><a href="mailto:sloisel@gmail.com" style="font-size: 13px;">sloisel@gmail.com</a><span style="font-size: 13px;">> wr=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Hi Barry,</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Thanks for your email. Thanks for pointing out my PetSC sillyness. I thin=</span><br style="font-size: 13px;"><span style="font-size: 13px;">k what happened is I played with matrices as well as with preconditioners s=</span><br style="font-size: 13px;"><span style="font-size: 13px;">o I initially implemented it as MATSHELL and at the end wrapped it in a PCS=</span><br style="font-size: 13px;"><span style="font-size: 13px;">HELL. :)</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> On Mon, Jan 23, 2017 at 2:26 AM, Barry Smith <</span><a href="mailto:bsmith@mcs.anl.gov" style="font-size: 13px;">bsmith@mcs.anl.gov</a><span style="font-size: 13px;">> wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> I've started to look at this; it is a little weird using MATSHELL within =</span><br style="font-size: 13px;"><span style="font-size: 13px;">your PCSHELL, seems unnecessary. Not necessarily wrong, just strange. I'll =</span><br style="font-size: 13px;"><span style="font-size: 13px;">continue to try to understand it.</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Barry</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > On Jan 20, 2017, at 5:48 PM, S=C3=A9bastien Loisel <</span><a href="mailto:sloisel@gmail.com" style="font-size: 13px;">sloisel@gmail.com</a><span style="font-size: 13px;">> =</span><br style="font-size: 13px;"><span style="font-size: 13px;">wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > OK I'm attaching the prototype.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > It won't be 100% plug-and-play for you because in principle T1 and T2 a=</span><br style="font-size: 13px;"><span style="font-size: 13px;">re built on top of "sub-preconditioners" (AMG for T1 and BILU for T2) and j=</span><br style="font-size: 13px;"><span style="font-size: 13px;">udging from the PetSC architecture elsewhere I must assume you're going to =</span><br style="font-size: 13px;"><span style="font-size: 13px;">want to expose those in a rational way. At present, I've hard-coded some su=</span><br style="font-size: 13px;"><span style="font-size: 13px;">b-preconditioners, and in particular for T2 I had to resort to PILUT becaus=</span><br style="font-size: 13px;"><span style="font-size: 13px;">e I didn't have a BILU(0) handy.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > Also, I broke the PetSC law and my functions don't return integers, so =</span><br style="font-size: 13px;"><span style="font-size: 13px;">I also assume you're going to want to tweak that... Sorry!</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > On Fri, Jan 20, 2017 at 11:29 PM, Barry Smith <</span><a href="mailto:bsmith@mcs.anl.gov" style="font-size: 13px;">bsmith@mcs.anl.gov</a><span style="font-size: 13px;">> wrot=</span><br style="font-size: 13px;"><span style="font-size: 13px;">e:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > Sure, email the PCSHELL version, best case scenario I only need chan=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ge out the PCSHELL and it takes me 5 minutes :-)</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > On Jan 20, 2017, at 5:07 PM, S=C3=A9bastien Loisel <</span><a href="mailto:sloisel@gmail.com" style="font-size: 13px;">sloisel@gmail.com</a><span style="font-size: 13px;">=</span><br style="font-size: 13px;"><span style="font-size: 13px;">> wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Hi all,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Thanks for your emails. I'm willing to help in whatever way. We have =</span><br style="font-size: 13px;"><span style="font-size: 13px;">a "PCSHELL" prototype we can provide on request, although PetSC experts can=</span><br style="font-size: 13px;"><span style="font-size: 13px;">no doubt do a better job than I did.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Thanks,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > On Fri, Jan 20, 2017 at 9:50 PM, Robert Annewandter <robert.annewandt=</span><br style="font-size: 13px;"><a href="mailto:er@opengosim.com" style="font-size: 13px;">er@opengosim.com</a><span style="font-size: 13px;">> wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Indeed that would be very helpful! We're very happy to support to tes=</span><br style="font-size: 13px;"><span style="font-size: 13px;">t things out, provide feedback etc</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Many thanks!</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Robert</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > On 20/01/17 21:22, Hammond, Glenn E wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >> That sounds great. I do know that there is also much interest state=</span><br style="font-size: 13px;"><span style="font-size: 13px;">-side in CPR preconditioning within PFLOTRAN. I have a Sandia colleague in=</span><br style="font-size: 13px;"><span style="font-size: 13px;">Carlsbad, NM who has been asking about it. I am sure that Sebastien and/o=</span><br style="font-size: 13px;"><span style="font-size: 13px;">r Robert will help out in any way possible.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >> Thanks,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >> Glenn</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> -----Original Message-----</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> From: Barry Smith [</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> </span><a href="mailto:bsmith@mcs.anl.gov" style="font-size: 13px;">mailto:bsmith@mcs.anl.gov</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> ]</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Sent: Friday, January 20, 2017 12:58 PM</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> To: Hammond, Glenn E</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> <</span><a href="mailto:gehammo@sandia.gov" style="font-size: 13px;">gehammo@sandia.gov</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Cc: S=C3=A9bastien Loisel</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> <</span><a href="mailto:sloisel@gmail.com" style="font-size: 13px;">sloisel@gmail.com</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> ; Robert Annewandter</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> <</span><a href="mailto:robert.annewandter@opengosim.com" style="font-size: 13px;">robert.annewandter@opengosim.com</a><span style="font-size: 13px;">>; Jed Brown <</span><a href="mailto:jed@jedbrown.org" style="font-size: 13px;">jed@jedbrown.org</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> ;</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Paolo Orsini</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> <</span><a href="mailto:paolo.orsini@opengosim.com" style="font-size: 13px;">paolo.orsini@opengosim.com</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> ; Matthew Knepley</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> <</span><a href="mailto:knepley@gmail.com" style="font-size: 13px;">knepley@gmail.com</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Subject: Re: [EXTERNAL] CPR preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Glenn,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Sorry about the delay in processing this, too much going on ...</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> I think the best thing is for us (the PETSc developers) to impl=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ement a CPR</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> preconditioner "directly" as its own PC and then have you guys try =</span><br style="font-size: 13px;"><span style="font-size: 13px;">it out. I am</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> planning to do this.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Barry</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> On Jan 20, 2017, at 2:50 PM, Hammond, Glenn E <</span><a href="mailto:gehammo@sandia.gov" style="font-size: 13px;">gehammo@sandia.gov</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> wrote:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Barry, Jed or Matt,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Do you have any suggestions for how to address the limitations of</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> PetscFieldSplit() discussed below. Will they need to manipulate th=</span><br style="font-size: 13px;"><span style="font-size: 13px;">e matrices</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> manually?</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Thanks,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Glenn</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> From: S=C3=A9bastien Loisel [</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> </span><a href="mailto:sloisel@gmail.com" style="font-size: 13px;">mailto:sloisel@gmail.com</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> ]</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Sent: Wednesday, January 11, 2017 3:33 AM</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> To: Barry Smith</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> <</span><a href="mailto:bsmith@mcs.anl.gov" style="font-size: 13px;">bsmith@mcs.anl.gov</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> ; Robert Annewandter</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> <</span><a href="mailto:robert.annewandter@opengosim.com" style="font-size: 13px;">robert.annewandter@opengosim.com</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> ; Hammond, Glenn E</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> <</span><a href="mailto:gehammo@sandia.gov" style="font-size: 13px;">gehammo@sandia.gov</a><span style="font-size: 13px;">>; Jed Brown <</span><a href="mailto:jed@jedbrown.org" style="font-size: 13px;">jed@jedbrown.org</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> ; Paolo Orsini</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> <</span><a href="mailto:paolo.orsini@opengosim.com" style="font-size: 13px;">paolo.orsini@opengosim.com</a><span style="font-size: 13px;">></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Subject: [EXTERNAL] CPR preconditioning</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Dear Friends,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Paolo has asked me to write this email to clarify issues surroundi=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ng</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> the CPR preconditioner that is widely used in multiphase flow. I k=</span><br style="font-size: 13px;"><span style="font-size: 13px;">now</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Barry from a long time ago but we only met once when I was a PhD</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> student so I would be shocked if he remembered me. :)</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> I'm a math assistant professor and one of my areas of specializati=</span><br style="font-size: 13px;"><span style="font-size: 13px;">on is linear</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> algebra and preconditioning.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> The main issue that is useful to clarify is the following. There w=</span><br style="font-size: 13px;"><span style="font-size: 13px;">as a proposal</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> to use PetSC's PETSCFIELDSPLIT in conjunction with PCCOMPOSITE in o=</span><br style="font-size: 13px;"><span style="font-size: 13px;">rder to</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> implement CPR preconditioning. Although this is morally the right i=</span><br style="font-size: 13px;"><span style="font-size: 13px;">dea, this</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> seems to be currently impossible because PETSCFIELDSPLIT lacks the</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> capabilities it would need to implement the T1 preconditioner. This=</span><br style="font-size: 13px;"><span style="font-size: 13px;">is due to</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> limitations in the API exposed by PETSCFIELDSPLIT (and no doubt lim=</span><br style="font-size: 13px;"><span style="font-size: 13px;">itations</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> in the underlying implementation).</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> In order to be as clear as possible, please allow me to describe</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> unambiguously the first of the two parts of the CPR preconditioner =</span><br style="font-size: 13px;"><span style="font-size: 13px;">using</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> some MATLAB snippets. Let I denote the N by N identity, and Z the N=</span><br style="font-size: 13px;"><span style="font-size: 13px;">by N</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> zero matrix. Put WT =3D [I I I] and C =3D [Z;Z;I]. The pressure mat=</span><br style="font-size: 13px;"><span style="font-size: 13px;">rix is Ap =3D</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> WT*A*C, and the T1 preconditioner is C*inv(Ap)*WT, where inv(Ap) is=</span><br style="font-size: 13px;"><span style="font-size: 13px;">to be</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> implemented with AMG.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> This T1 preconditioner is the one that would have to be implemente=</span><br style="font-size: 13px;"><span style="font-size: 13px;">d by</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> PETSCFIELDSPLIT. The following limitations in PETSCFIELDSPLIT preve=</span><br style="font-size: 13px;"><span style="font-size: 13px;">nts one</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> to implement T1:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> =E2=80=A2 One must select the "pressure" by specifying an IS ob=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ject to</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> PCFieldSplitSetIS(). However, since WT =3D [I I I], the pressure is=</span><br style="font-size: 13px;"><span style="font-size: 13px;">obtained by</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> summing over the three fields. As far as I can tell, an IS object d=</span><br style="font-size: 13px;"><span style="font-size: 13px;">oes not allow</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> one to sum over several entries to obtain the pressure field.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> =E2=80=A2 The pressure matrix is Ap =3D WT*A*C; note that the m=</span><br style="font-size: 13px;"><span style="font-size: 13px;">atrix WT on</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> the left is different from the matrix C on the right. However, PCFI=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ELDSPLIT</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> has no notion of a "left-IS" vs "right-IS"; morally, only diagonal =</span><br style="font-size: 13px;"><span style="font-size: 13px;">blocks of A can</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> be used by PCFIELDSPLIT.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> =E2=80=A2 PCFIELDSPLIT offers a range of hard-coded block struc=</span><br style="font-size: 13px;"><span style="font-size: 13px;">tures for the</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> final preconditioner, but the choice T1 =3D C*inv(Ap)*WT is not one=</span><br style="font-size: 13px;"><span style="font-size: 13px;">of these</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> choices. Indeed, this first stage CPR preconditioner T1 is *singula=</span><br style="font-size: 13px;"><span style="font-size: 13px;">r*, but there</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> is no obvious way for PCFIELDSPLIT to produce a singular preconditi=</span><br style="font-size: 13px;"><span style="font-size: 13px;">oner.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Note that the documentation for PETSCFIELDSPLIT says that "The</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> Constrained Pressure Preconditioner (CPR) does not appear to be cur=</span><br style="font-size: 13px;"><span style="font-size: 13px;">rently</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> implementable directly with PCFIELDSPLIT". Unless there are very si=</span><br style="font-size: 13px;"><span style="font-size: 13px;">gnificant</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> capabilities that are not documented, I don't see how CPR can be</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> implemented with PETSCFIELDSPLIT.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Elsewhere, someone proposed putting the two preconditioners T1 and=</span><br style="font-size: 13px;"><span style="font-size: 13px;">T2</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> on each side of A, e.g. T1*A*T2. That is a very bad idea because T=</span><br style="font-size: 13px;"><span style="font-size: 13px;">1 is</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> singular and hence T1*A*T2 is also singular. The correct CPR</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> preconditioner is nonsingular despite the fact that T1 is singular,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> and MCPR is given by the formula MCPR =3D T2*(I-A*T1)+T1, where T2=</span><br style="font-size: 13px;"><span style="font-size: 13px;">=3D</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> BILU(0) of A. (There is a proof, due to Felix Kwok, that BILU(0) w=</span><br style="font-size: 13px;"><span style="font-size: 13px;">orks</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> even though ILU(0) craps out on vanishing diagonal entries.)</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> I'm also attaching a sample MATLAB code that runs the CPR precondi=</span><br style="font-size: 13px;"><span style="font-size: 13px;">tioner</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> on some fabricated random matrix A. I emphasize that this is not a =</span><br style="font-size: 13px;"><span style="font-size: 13px;">realistic</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> matrix, but it still illustrates that the algorithm works, and that=</span><br style="font-size: 13px;"><span style="font-size: 13px;">MCPR is better</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> than T2 alone. Note that T1 cannot be used alone since it is singul=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ar. Further</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> gains are expected when the Robert's realistic code with correct ph=</span><br style="font-size: 13px;"><span style="font-size: 13px;">ysics will</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>> come online.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> <image003.jpg></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> I hope this clarifies some things.</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> S=C3=A9bastien Loisel</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Assistant Professor</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> Department of Mathematics, Heriot-Watt University Riccarton, EH14 =</span><br style="font-size: 13px;"><span style="font-size: 13px;">4AS,</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> United Kingdom</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> web:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> </span><a href="http://www.ma.hw.ac.uk/~loisel/" style="font-size: 13px;">http://www.ma.hw.ac.uk/~loisel/</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> email: S.Loisel at</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> </span><a href="http://hw.ac.uk/" style="font-size: 13px;">hw.ac.uk</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> phone:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> +44 131 451 3234</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> fax:</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > >>>> +44 131 451 3249</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > --</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > S=C3=A9bastien Loisel</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Assistant Professor</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Department of Mathematics, Heriot-Watt University</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > Riccarton, EH14 4AS, United Kingdom</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > web: </span><a href="http://www.ma.hw.ac.uk/~loisel/" style="font-size: 13px;">http://www.ma.hw.ac.uk/~loisel/</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > > email: S.Loisel at </span><a href="http://hw.ac.uk/" style="font-size: 13px;">hw.ac.uk</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > > phone: +44 131 451 3234</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > > fax: +44 131 451 3249</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > --</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > S=C3=A9bastien Loisel</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > Assistant Professor</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > Department of Mathematics, Heriot-Watt University</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > Riccarton, EH14 4AS, United Kingdom</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > web: </span><a href="http://www.ma.hw.ac.uk/~loisel/" style="font-size: 13px;">http://www.ma.hw.ac.uk/~loisel/</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > email: S.Loisel at </span><a href="http://hw.ac.uk/" style="font-size: 13px;">hw.ac.uk</a><br style="font-size: 13px;"><span style="font-size: 13px;">> > phone: +44 131 451 3234</span><br style="font-size: 13px;"><span style="font-size: 13px;">> > fax: +44 131 451 3249</span><br style="font-size: 13px;"><span style="font-size: 13px;">> ></span><br style="font-size: 13px;"><span style="font-size: 13px;">> > <petsc.zip></span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> --=20</span><br style="font-size: 13px;"><span style="font-size: 13px;">> S=C3=A9bastien Loisel</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Assistant Professor</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Department of Mathematics, Heriot-Watt University</span><br style="font-size: 13px;"><span style="font-size: 13px;">> Riccarton, EH14 4AS, United Kingdom</span><br style="font-size: 13px;"><span style="font-size: 13px;">> web: </span><a href="http://www.ma.hw.ac.uk/~loisel/" style="font-size: 13px;">http://www.ma.hw.ac.uk/~loisel/</a><br style="font-size: 13px;"><span style="font-size: 13px;">> email: S.Loisel at </span><a href="http://hw.ac.uk/" style="font-size: 13px;">hw.ac.uk</a><br style="font-size: 13px;"><span style="font-size: 13px;">> phone: +44 131 451 3234</span><br style="font-size: 13px;"><span style="font-size: 13px;">> fax: +44 131 451 3249</span><br style="font-size: 13px;"><span style="font-size: 13px;">>=20</span><br style="font-size: 13px;"><br style="font-size: 13px;"><br style="font-size: 13px;"><span style="font-size: 13px;">--====-=-=--</span><br style="font-size: 13px;"><span style="font-size: 13px;">]</span></div></body></html>