PCBlockDiagonal.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 "PCBlockDiagonal.hpp"
00033 #include "Exception.hpp"
00034 
00035 PCBlockDiagonal::PCBlockDiagonal(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     PCBlockDiagonalCreate(rKspObject);
00045     PCBlockDiagonalSetUp();
00046 }
00047 
00048 PCBlockDiagonal::~PCBlockDiagonal()
00049 {
00050 #ifdef TRACE_KSP
00051     if (PetscTools::AmMaster())
00052     {
00053         std::cout << " -- Block diagonal 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 mGatherTime: " << mPCContext.mGatherTime << std::endl;
00058     }
00059 #endif
00060 
00061     MatDestroy(mPCContext.A11_matrix_subblock);
00062     MatDestroy(mPCContext.A22_matrix_subblock);
00063 
00064     PCDestroy(mPCContext.PC_amg_A11);
00065     PCDestroy(mPCContext.PC_amg_A22);
00066 
00067     VecDestroy(mPCContext.x1_subvector);
00068     VecDestroy(mPCContext.y1_subvector);
00069 
00070     VecDestroy(mPCContext.x2_subvector);
00071     VecDestroy(mPCContext.y2_subvector);
00072     
00073     VecScatterDestroy(mPCContext.A11_scatter_ctx);
00074     VecScatterDestroy(mPCContext.A22_scatter_ctx);    
00075 }
00076 
00077 void PCBlockDiagonal::PCBlockDiagonalCreate(KSP& rKspObject)
00078 {
00079     KSPGetPC(rKspObject, &mPetscPCObject);
00080 
00081     Mat system_matrix, dummy;
00082     MatStructure flag;
00083     KSPGetOperators(rKspObject, &system_matrix, &dummy, &flag);
00084 
00085     PetscInt num_rows, num_columns;
00086     MatGetSize(system_matrix, &num_rows, &num_columns);
00087     assert(num_rows==num_columns);
00088 
00089     PetscInt num_local_rows, num_local_columns;
00090     MatGetLocalSize(system_matrix, &num_local_rows, &num_local_columns);
00091 
00092     // odd number of rows: impossible in Bidomain.
00093     // odd number of local rows: impossible if V_m and phi_e for each node are stored in the same processor.
00094     if ((num_rows%2 != 0) || (num_local_rows%2 != 0))
00095     {
00096         TERMINATE("Wrong matrix parallel layout detected in PCLDUFactorisation.");
00097     }
00098 
00099     // Allocate memory     
00100     unsigned subvector_num_rows = num_rows/2;    
00101     unsigned subvector_local_rows = num_local_rows/2;
00102     mPCContext.x1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00103     mPCContext.x2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00104     mPCContext.y1_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00105     mPCContext.y2_subvector = PetscTools::CreateVec(subvector_num_rows, subvector_local_rows);
00106 
00107     // Create scatter contexts
00108     {
00109         // Needed by SetupInterleavedVectorScatterGather in order to find out parallel layout.
00110         Vec dummy_vec = PetscTools::CreateVec(num_rows, num_local_rows);
00111 
00112         PetscVecTools::SetupInterleavedVectorScatterGather(dummy_vec, mPCContext.A11_scatter_ctx, mPCContext.A22_scatter_ctx);
00113 
00114         VecDestroy(dummy_vec);        
00115     }
00116             
00117     // Get matrix sublock A11        
00118     {
00119         // Work out local row range for subblock A11 (same as x1 or y1) 
00120         PetscInt low, high, global_size;
00121         VecGetOwnershipRange(mPCContext.x1_subvector, &low, &high);
00122         VecGetSize(mPCContext.x1_subvector, &global_size);        
00123         assert(global_size == num_rows/2);       
00124         
00125         IS A11_local_rows;
00126         IS A11_columns;
00127         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low, 2, &A11_local_rows);
00128         ISCreateStride(PETSC_COMM_WORLD, global_size, 0, 2, &A11_columns);
00129     
00130 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00131         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns,
00132             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00133 #else
00134         MatGetSubMatrix(system_matrix, A11_local_rows, A11_columns, PETSC_DECIDE, 
00135             MAT_INITIAL_MATRIX, &mPCContext.A11_matrix_subblock);
00136 #endif
00137 
00138     
00139         ISDestroy(A11_local_rows);
00140         ISDestroy(A11_columns);
00141     }
00142 
00143     // Get matrix sublock A22
00144     {
00145         // Work out local row range for subblock A22 (same as x2 or y2) 
00146         PetscInt low, high, global_size;
00147         VecGetOwnershipRange(mPCContext.x2_subvector, &low, &high);
00148         VecGetSize(mPCContext.x2_subvector, &global_size);        
00149         assert(global_size == num_rows/2);       
00150         
00151         IS A22_local_rows;
00152         IS A22_columns;
00153         ISCreateStride(PETSC_COMM_WORLD, high-low, 2*low+1, 2, &A22_local_rows);
00154         ISCreateStride(PETSC_COMM_WORLD, global_size, 1, 2, &A22_columns);
00155     
00156 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00157         MatGetSubMatrix(system_matrix, A22_local_rows, A22_columns,
00158             MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
00159 #else
00160         MatGetSubMatrix(system_matrix, A22_local_rows, A22_columns, PETSC_DECIDE, 
00161             MAT_INITIAL_MATRIX, &mPCContext.A22_matrix_subblock);
00162 #endif
00163     
00164         ISDestroy(A22_local_rows);
00165         ISDestroy(A22_columns);
00166     }
00167 
00168     // Register call-back function and its context
00169     PCSetType(mPetscPCObject, PCSHELL);
00170 #if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
00171     PCShellSetApply(mPetscPCObject, PCBlockDiagonalApply, (void*) &mPCContext);
00172 #else
00173     // Register PC context so it gets passed to PCBlockDiagonalApply
00174     PCShellSetContext(mPetscPCObject, &mPCContext);
00175     // Register call-back function
00176     PCShellSetApply(mPetscPCObject, PCBlockDiagonalApply);
00177 #endif
00178 
00179 }
00180 
00181 void PCBlockDiagonal::PCBlockDiagonalSetUp()
00182 {
00183     // These options will get read by PCSetFromOptions
00184 //     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00185 //     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00186 //     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00187 
00188     // Set up amg preconditioner for block A11
00189     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A11));
00190 
00191     //    PCSetType(mPCContext.PC_amg_A11, PCBJACOBI);
00192 
00193     PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
00194     //PCHYPRESetType(mPCContext.PC_amg_A11, "euclid");
00195     PetscOptionsSetValue("-pc_hypre_type", "euclid");
00196     PetscOptionsSetValue("-pc_hypre_euclid_levels", "0");
00197 
00198 
00199     //PCSetType(mPCContext.PC_amg_A11, PCHYPRE);
00200     //PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00201     //PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00202     //PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00203     //PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type", "HMIS");
00204 
00205     PCSetOperators(mPCContext.PC_amg_A11, mPCContext.A11_matrix_subblock, mPCContext.A11_matrix_subblock, SAME_PRECONDITIONER);
00206     PCSetFromOptions(mPCContext.PC_amg_A11);
00207     PCSetUp(mPCContext.PC_amg_A11);
00208 
00209     // Set up amg preconditioner for block A22
00210     PCCreate(PETSC_COMM_WORLD, &(mPCContext.PC_amg_A22));
00211 
00212     /* Full AMG in the block */
00213     PCSetType(mPCContext.PC_amg_A22, PCHYPRE);
00214     //PCHYPRESetType(mPCContext.PC_amg_A22, "boomeramg");
00215     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00216     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00217     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00218     PetscOptionsSetValue("-pc_hypre_boomeramg_coarsen_type", "HMIS");
00219     //    PetscOptionsSetValue("-pc_hypre_boomeramg_relax_type_all","Jacobi");
00220     //PetscOptionsSetValue("-pc_hypre_boomeramg_max_levels","10");
00221     //PetscOptionsSetValue("-pc_hypre_boomeramg_agg_nl", "1");
00222     //    PetscOptionsSetValue("-pc_hypre_boomeramg_print_statistics","");
00223     //    PetscOptionsSetValue("-pc_hypre_boomeramg_interp_type","standard");
00224     //    PetscOptionsSetValue("-pc_hypre_boomeramg_interp_type","ext+i");
00225 
00226     //    PCHYPRESetType(mPCContext.PC_amg_A22, "euclid");
00227 
00228 
00229 
00230     /* Block Jacobi with AMG at each block */
00231     //     PCSetType(mPCContext.PC_amg_A22, PCBJACOBI);
00232     
00233     //     PetscOptionsSetValue("-sub_pc_type", "hypre");
00234     
00235     //     PetscOptionsSetValue("-sub_pc_hypre_type", "boomeramg");
00236     //     PetscOptionsSetValue("-sub_pc_hypre_boomeramg_max_iter", "1");
00237     //     PetscOptionsSetValue("-sub_pc_hypre_boomeramg_strong_threshold", "0.0");
00238     
00239     //     PetscOptionsSetValue("-pc_hypre_type", "boomeramg");
00240     //     PetscOptionsSetValue("-pc_hypre_boomeramg_max_iter", "1");
00241     //     PetscOptionsSetValue("-pc_hypre_boomeramg_strong_threshold", "0.0");
00242 
00243     /* Additive Schwarz with AMG at each block */
00244 //     PCSetType(mPCContext.PC_amg_A22, PCASM);
00245     
00246 //     PetscOptionsSetValue("-pc_asm_type", "basic");
00247 //     PetscOptionsSetValue("-pc_asm_overlap", "1");
00248 
00249 //     PetscOptionsSetValue("-sub_ksp_type", "preonly");    
00250 
00251 //     PetscOptionsSetValue("-sub_pc_type", "hypre");
00252     
00253 //     PetscOptionsSetValue("-sub_pc_hypre_type", "boomeramg");
00254 //     PetscOptionsSetValue("-sub_pc_hypre_boomeramg_max_iter", "1");
00255 //     PetscOptionsSetValue("-sub_pc_hypre_boomeramg_strong_threshold", "0.0");
00256  
00257 
00258     PCSetOperators(mPCContext.PC_amg_A22, mPCContext.A22_matrix_subblock, mPCContext.A22_matrix_subblock, SAME_PRECONDITIONER);
00259     PCSetFromOptions(mPCContext.PC_amg_A22);
00260     PCSetUp(mPCContext.PC_amg_A22);
00261 }
00262 
00263 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 1) //PETSc 3.1
00264 PetscErrorCode PCBlockDiagonalApply(PC pc_object, Vec x, Vec y)
00265 {
00266   void* pc_context;
00267 
00268   PCShellGetContext(pc_object, &pc_context);   
00269 #else
00270 PetscErrorCode PCBlockDiagonalApply(void* pc_context, Vec x, Vec y)
00271 {
00272 #endif
00273 
00274     // Cast the context pointer to PCBlockDiagonalContext
00275     PCBlockDiagonal::PCBlockDiagonalContext* block_diag_context = (PCBlockDiagonal::PCBlockDiagonalContext*) pc_context;
00276     assert(block_diag_context!=NULL);
00277 
00278     /*
00279      * Scatter x = [x1 x2]'
00280      */
00281 #ifdef TRACE_KSP
00282     double init_time = MPI_Wtime();
00283 #endif
00284 
00285     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);
00286 
00287 #ifdef TRACE_KSP
00288     block_diag_context->mScatterTime += MPI_Wtime() - init_time;
00289 #endif
00290 
00291     /*
00292      *  y1 = AMG(A11)*x1
00293      *  y2 = AMG(A22)*x2
00294      */
00295 #ifdef TRACE_KSP
00296     init_time = MPI_Wtime();
00297 #endif
00298     PCApply(block_diag_context->PC_amg_A11, block_diag_context->x1_subvector, block_diag_context->y1_subvector);
00299 #ifdef TRACE_KSP
00300     block_diag_context->mA1PreconditionerTime += MPI_Wtime() - init_time;
00301 #endif
00302 
00303 #ifdef TRACE_KSP
00304     init_time = MPI_Wtime();
00305 #endif
00306     PCApply(block_diag_context->PC_amg_A22, block_diag_context->x2_subvector, block_diag_context->y2_subvector);
00307 #ifdef TRACE_KSP
00308     block_diag_context->mA2PreconditionerTime += MPI_Wtime() - init_time;
00309 #endif
00310 
00311     /*
00312      * Gather y = [y1 y2]'
00313      */
00314 //PETSc-3.x.x or PETSc-2.3.3
00315 #ifdef TRACE_KSP
00316     init_time = MPI_Wtime();
00317 #endif
00318 
00319     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);
00320 
00321 #ifdef TRACE_KSP
00322     block_diag_context->mGatherTime += MPI_Wtime() - init_time;
00323 #endif
00324     return 0;
00325 }

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