PCLDUFactorisation.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 #include <iostream>
00030 
00031 #include "PetscVecTools.hpp" // Includes Ublas so must come first
00032 #include "PCLDUFactorisation.hpp"
00033 #include "Exception.hpp"
00034 
00035 PCLDUFactorisation::PCLDUFactorisation(KSP& rKspObject)
00036 {
00037 #ifdef TRACE_KSP
00038     mPCContext.mScatterTime = 0.0;
00039     mPCContext.mA1PreconditionerTime = 0.0;;
00040     mPCContext.mA2PreconditionerTime = 0.0;;
00041     mPCContext.mGatherTime = 0.0;;
00042 #endif
00043 
00044     PCLDUFactorisationCreate(rKspObject);
00045     PCLDUFactorisationSetUp();
00046 }
00047 
00048 PCLDUFactorisation::~PCLDUFactorisation()
00049 {
00050 #ifdef TRACE_KSP
00051     if (PetscTools::AmMaster())
00052     {
00053         std::cout << " -- LDU factorisation preconditioner profile information: " << std::endl;
00054         std::cout << "\t mScatterTime: " << mPCContext.mScatterTime << std::endl;
00055         std::cout << "\t mA1PreconditionerTime: " << mPCContext.mA1PreconditionerTime << std::endl;
00056         std::cout << "\t mA2PreconditionerTime: " << mPCContext.mA2PreconditionerTime << std::endl;
00057         std::cout << "\t mExtraLAOperations: " << mPCContext.mExtraLAOperations << std::endl;
00058         std::cout << "\t mGatherTime: " << mPCContext.mGatherTime << std::endl;
00059     }
00060 #endif
00061 
00062     MatDestroy(mPCContext.A11_matrix_subblock);
00063     MatDestroy(mPCContext.A22_matrix_subblock);
00064     MatDestroy(mPCContext.B_matrix_subblock);
00065 
00066     PCDestroy(mPCContext.PC_amg_A11);
00067     PCDestroy(mPCContext.PC_amg_A22);
00068 
00069     VecDestroy(mPCContext.x1_subvector);
00070     VecDestroy(mPCContext.y1_subvector);
00071     VecDestroy(mPCContext.x2_subvector);
00072     VecDestroy(mPCContext.y2_subvector);
00073     VecDestroy(mPCContext.z);
00074     VecDestroy(mPCContext.temp);
00075     
00076     VecScatterDestroy(mPCContext.A11_scatter_ctx);
00077     VecScatterDestroy(mPCContext.A22_scatter_ctx);        
00078 }
00079 
00080 void PCLDUFactorisation::PCLDUFactorisationCreate(KSP& rKspObject)
00081 {
00082     KSPGetPC(rKspObject, &mPetscPCObject);
00083 
00084     Mat system_matrix, dummy;
00085     MatStructure flag;
00086     KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
00087 
00088     PetscInt num_rows, num_columns;
00089     MatGetSize(system_matrix, &num_rows, &num_columns);
00090 
00091     PetscInt num_local_rows, num_local_columns;
00092     MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
00093 
00094     // odd number of rows: impossible in Bidomain.
00095     // odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
00096     if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
00097     {
00098         TERMINATE("Wrong matrix parallel layout detected in PCLDUFactorisation.");
00099     }
00100 
00101     // Allocate memory     
00102     unsigned subvector_num_rows = num_rows/2;    
00103     unsigned subvector_local_rows = num_local_rows/2;
00104     mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00105     mPCContext.x2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00106     mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00107     mPCContext.y2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00108     mPCContext.z = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00109     mPCContext.temp = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00110 
00111     // Create scatter contexts
00112     {
00113         // Needed by VecScatterCreate in order to find out parallel layout.
00114         Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
00115 
00116         PetscVecTools::SetupInterleavedVectorScatterGather(dummy_vec, mPCContext.A11_scatter_ctx, mPCContext.A22_scatter_ctx);
00117 
00118         VecDestroy(dummy_vec);
00119     }
00120 
00121 
00122     // Get matrix sublock A11        
00123     {
00124         // Work out local row range for subblock A11 (same as x1 or y1) 
00125         PetscInt low, high, global_size;
00126         VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
00127         VecGetSize(mPCContext.x1_subvector, &global_size);        
00128         assert(global_size == num_rows/2);       
00129         
00130         IS A11_local_rows;
00131         IS A11_columns;
00132         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows);
00133         ISCreateStride(PETSC_COMM_WORLD, global_size, 0, 2, &A11_columns);
00134     
00135 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00136         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns,
00137             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00138 #else
00139         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns, PETSC_DECIDE,
00140             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00141 #endif    
00142         ISDestroy(A11_local_rows);
00143         ISDestroy(A11_columns);
00144     }
00145 
00146     // Get matrix sublock A22
00147     {
00148         // Work out local row range for subblock A22 (same as x2 or y2) 
00149         PetscInt low, high, global_size;
00150         VecGetOwnershipRange(mPCContext.x2_subvector, &low, &high);
00151         VecGetSize(mPCContext.x2_subvector, &global_size);        
00152         assert(global_size == num_rows/2);       
00153         
00154         IS A22_local_rows;
00155         IS A22_columns;
00156         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low+1, 2, &A22_local_rows);
00157         ISCreateStride(PETSC_COMM_WORLD, global_size, 1, 2, &A22_columns);
00158 
00159 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00160         MatGetSubMatrix(system_matrix, A22_local_rows, A22_columns,
00161             MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
00162 #else
00163         MatGetSubMatrix(system_matrix, A22_local_rows, A22_columns, PETSC_DECIDE, 
00164             MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
00165 #endif
00166     
00167         ISDestroy(A22_local_rows);
00168         ISDestroy(A22_columns);
00169     }
00170 
00171     // Get matrix sublock B (the upper triangular one)
00172     {
00173         // Work out local row range for subblock B (same as A11, x1 or y1) 
00174         PetscInt low, high, global_size;
00175         VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
00176         VecGetSize(mPCContext.x1_subvector, &global_size);        
00177         assert(global_size == num_rows/2);        
00178         
00179         IS B_local_rows;
00180         IS B_columns;
00181         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &B_local_rows);
00182         ISCreateStride(PETSC_COMM_WORLD, global_size, 1, 2, &B_columns);
00183     
00184 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00185         MatGetSubMatrix(system_matrix, B_local_rows, B_columns,
00186             MAT_INITIAL_MATRIX, &mPCContext.B_matrix_subblock);        
00187 #else
00188         MatGetSubMatrix(system_matrix, B_local_rows, B_columns, PETSC_DECIDE, 
00189             MAT_INITIAL_MATRIX, &mPCContext.B_matrix_subblock);
00190 #endif
00191     
00192         ISDestroy(B_local_rows);
00193         ISDestroy(B_columns);
00194     }
00195     
00196     /*
00197      * Experimental (#1082): in PP removing the mass matrix from the A22 block seems to work better.
00198      *                       This is equivalent to do A22 = A22 + B in this implementation. 
00199      */
00200 // #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00201 //     PetscScalar petsc_one = 1.0;
00202 //     MatAXPY(&petsc_one, mPCContext.B_matrix_subblock, mPCContext.A22_matrix_subblock, DIFFERENT_NONZERO_PATTERN);
00203 // #else
00204 //     MatAXPY(mPCContext.A22_matrix_subblock, 1.0, mPCContext.B_matrix_subblock, DIFFERENT_NONZERO_PATTERN);
00205 // #endif    
00206     
00207 //     // Shift the block
00208 // #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00209 //     PetscScalar shift = -1e-8;
00210 //     MatShift(&shift, mPCContext.A22_matrix_subblock);
00211 // #else
00212 //     PetscScalar shift = -1e-8;
00213 //     MatShift(mPCContext.A22_matrix_subblock, shift);
00214 // #endif    
00215 
00216 
00217 
00218     PCSetType(mPetscPCObject, PCSHELL);
00219 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00220     // Register PC context and call-back function
00221     PCShellSetApply(mPetscPCObject, PCLDUFactorisationApply, (void*) &mPCContext);
00222 #else
00223     // Register PC context so it gets passed to PCBlockDiagonalApply
00224     PCShellSetContext(mPetscPCObject, &mPCContext);
00225     // Register call-back function
00226     PCShellSetApply(mPetscPCObject, PCLDUFactorisationApply);
00227 #endif
00228 
00229 }
00230 
00231 void PCLDUFactorisation::PCLDUFactorisationSetUp()
00232 {
00233     // These options will get read by PCSetFromOptions
00234 //     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00235 //     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00236 //     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00237 
00238     /*
00239      *  Set up preconditioner for block A11
00240      */
00241     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
00242     PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, SAME_PRECONDITIONER);
00243 
00244     // Choose between the two following blocks in order to approximate inv(A11) with one AMG cycle
00245     // or with an CG solve with high tolerance
00247 //    PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
00248 
00249     PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
00250     //PCHYPRESetType(mPCContext.PC_amg_A11, "euclid");
00251     PetscOptionsSetValue("-pc_hypre_type", "euclid");
00252     PetscOptionsSetValue("-pc_hypre_euclid_levels", "0");
00253     
00254 //     PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
00255 //     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00256 //     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00257 //     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00258 //     PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type", "HMIS");
00259 
00261 //    PCSetType(mPCContext.PC_amg_A11, PCKSP);
00262 //    KSP ksp1;
00263 //    PCKSPGetKSP(mPCContext.PC_amg_A11,&ksp1);
00264 //    KSPSetType(ksp1, KSPCG);
00265 //    KSPSetTolerances(ksp1, 0.1, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
00266 //
00267 //    PC prec1;
00268 //    KSPGetPC(ksp1, &prec1);
00269 //    PCSetType(prec1, PCBJACOBI);
00270 //    PCSetFromOptions(prec1);
00271 //    PCSetOperators(prec1, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, DIFFERENT_NONZERO_PATTERN);//   SAME_PRECONDITIONER);
00272 //    PCSetUp(prec1);
00273 //
00274 //    KSPSetFromOptions(ksp1);
00275 //    KSPSetUp(ksp1);
00277 
00278     PCSetFromOptions(mPCContext.PC_amg_A11);
00279     PCSetUp(mPCContext.PC_amg_A11);
00280 
00281 
00282     /*
00283      *  Set up amg preconditioner for block A22
00284      */
00285     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22));
00286     PCSetOperators(mPCContext.PC_amg_A22, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock, SAME_PRECONDITIONER);
00287 
00288     // Choose between the two following blocks in order to approximate inv(A11) with one AMG cycle
00289     // or with an CG solve with high tolerance
00291     PCSetType(mPCContext.PC_amg_A22, PCHYPRE);
00292     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00293     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00294     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00295     PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type", "HMIS");
00296     //    PetscOptionsSetValue("-pc_hypre_boomeramg_interp_type","ext+i");
00297 
00299 //    PCSetType(mPCContext.PC_amg_A22, PCKSP);
00300 //    KSP ksp2;
00301 //    PCKSPGetKSP(mPCContext.PC_amg_A22,&ksp2);
00302 //    KSPSetType(ksp2, KSPCG);
00303 //    KSPSetTolerances(ksp2, 0.1, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
00304 //
00305 //    PC prec2;
00306 //    KSPGetPC(ksp2, &prec2);
00307 //    PCSetType(prec2, PCBJACOBI);
00308 //    PCSetFromOptions(prec2);
00309 //    PCSetOperators(prec2, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock, DIFFERENT_NONZERO_PATTERN);//   SAME_PRECONDITIONER);
00310 //    PCSetUp(prec2);
00311 //
00312 //    KSPSetFromOptions(ksp2);
00313 //    KSPSetUp(ksp2);
00315 
00316     PCSetFromOptions(mPCContext.PC_amg_A22);
00317     PCSetUp(mPCContext.PC_amg_A22);
00318 }
00319 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00320 PetscErrorCode PCLDUFactorisationApply(PC pc_object, Vec x, Vec y)
00321 {
00322   void* pc_context;
00323 
00324   PCShellGetContext(pc_object, &pc_context);   
00325 #else
00326 PetscErrorCode PCLDUFactorisationApply(void* pc_context, Vec x, Vec y)
00327 {
00328 #endif
00330 
00331     // Cast the pointer to a PC context to our defined type
00332     PCLDUFactorisation::PCLDUFactorisationContext* block_diag_context = (PCLDUFactorisation::PCLDUFactorisationContext*) pc_context;
00333     assert(block_diag_context!=NULL);
00334 
00335     /*
00336      *  Split vector x into two. x = [x1 x2]'
00337      */
00338 #ifdef TRACE_KSP
00339     double init_time = MPI_Wtime();
00340 #endif
00341 
00342     PetscVecTools::DoInterleavedVecScatter(x, block_diag_context->A11_scatter_ctx, block_diag_context->x1_subvector, block_diag_context->A22_scatter_ctx, block_diag_context->x2_subvector);
00343 
00344 #ifdef TRACE_KSP
00345     block_diag_context->mScatterTime += MPI_Wtime() - init_time;
00346 #endif
00347 
00348     /*
00349      *  Apply preconditioner: [y1 y2]' = inv(P)[x1 x2]' 
00350      *
00351      *     z  = inv(A11)*x1
00352      *     y2 = inv(A22)*(x2 - B*z)
00353      *     y1 = z - inv(A11)(B*y2)
00354      */
00355 #ifdef TRACE_KSP
00356     init_time = MPI_Wtime();
00357 #endif
00358     //z  = inv(A11)*x1
00359     PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->z);
00360 #ifdef TRACE_KSP
00361     block_diag_context->mA1PreconditionerTime += MPI_Wtime() - init_time;
00362 #endif
00363 
00364     //y2 = inv(A22)*(x2 - B*z)
00365 #ifdef TRACE_KSP
00366     init_time = MPI_Wtime();
00367 #endif
00368     MatMult(block_diag_context->B_matrix_subblock,block_diag_context->z,block_diag_context->temp); //temp = B*z
00369     double minus_one = -1.0;
00370 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00371     VecAYPX(&minus_one, block_diag_context->x2_subvector, block_diag_context->temp); // temp <-- x2 - temp
00372 #else
00373     VecAYPX(block_diag_context->temp, minus_one, block_diag_context->x2_subvector); // temp <-- x2 - temp
00374 #endif
00375 #ifdef TRACE_KSP
00376     block_diag_context->mExtraLAOperations += MPI_Wtime() - init_time;
00377 #endif
00378 
00379 
00380 #ifdef TRACE_KSP
00381     init_time = MPI_Wtime();
00382 #endif
00383     PCApply(block_diag_context->PC_amg_A22, block_diag_context->temp, block_diag_context->y2_subvector); // y2 = inv(A22)*temp
00384 #ifdef TRACE_KSP
00385     block_diag_context->mA2PreconditionerTime += MPI_Wtime() - init_time;
00386 #endif
00387 
00388 
00389     //y1 = z - inv(A11)(B*y2)
00390 #ifdef TRACE_KSP
00391     init_time = MPI_Wtime();
00392 #endif
00393     MatMult(block_diag_context->B_matrix_subblock,block_diag_context->y2_subvector,block_diag_context->temp); //temp = B*y2
00394 #ifdef TRACE_KSP
00395     block_diag_context->mExtraLAOperations += MPI_Wtime() - init_time;
00396 #endif
00397 #ifdef TRACE_KSP
00398     init_time = MPI_Wtime();
00399 #endif
00400     PCApply(block_diag_context->PC_amg_A11, block_diag_context->temp, block_diag_context->y1_subvector); // y1 = inv(A11)*temp
00401 #ifdef TRACE_KSP
00402     block_diag_context->mA1PreconditionerTime += MPI_Wtime() - init_time;
00403 #endif
00404 
00405 #ifdef TRACE_KSP
00406     init_time = MPI_Wtime();
00407 #endif
00408 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00409     VecAYPX(&minus_one, block_diag_context->z, block_diag_context->y1_subvector); // y1 <-- z - y1
00410 #else
00411     VecAYPX(block_diag_context->y1_subvector, minus_one, block_diag_context->z); // y1 <-- z - y1
00412 #endif
00413 #ifdef TRACE_KSP
00414     block_diag_context->mExtraLAOperations += MPI_Wtime() - init_time;
00415 #endif
00416 
00417 
00418     /*
00419      *  Gather vectors y1 and y2. y = [y1 y2]'
00420      */
00421 #ifdef TRACE_KSP
00422     init_time = MPI_Wtime();
00423 #endif
00424 
00425     PetscVecTools::DoInterleavedVecGather(y, block_diag_context->A11_scatter_ctx, block_diag_context->y1_subvector, block_diag_context->A22_scatter_ctx, block_diag_context->y2_subvector);
00426 
00427 #ifdef TRACE_KSP
00428     block_diag_context->mGatherTime += MPI_Wtime() - init_time;
00429 #endif
00430 
00431     return 0;
00432 }

Generated on Tue May 31 14:31:42 2011 for Chaste by  doxygen 1.5.5