1 #include "bddc.h" 2 #include "bddcprivate.h" 3 #include <petscblaslapack.h> 4 5 #undef __FUNCT__ 6 #define __FUNCT__ "PCBDDCSetLevel" 7 PetscErrorCode PCBDDCSetLevel(PC pc,PetscInt level) 8 { 9 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 10 11 PetscFunctionBegin; 12 pcbddc->current_level=level; 13 PetscFunctionReturn(0); 14 } 15 16 #undef __FUNCT__ 17 #define __FUNCT__ "PCBDDCResetCustomization" 18 PetscErrorCode PCBDDCResetCustomization(PC pc) 19 { 20 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 21 PetscInt i; 22 PetscErrorCode ierr; 23 24 PetscFunctionBegin; 25 ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr); 26 ierr = ISDestroy(&pcbddc->user_primal_vertices);CHKERRQ(ierr); 27 ierr = MatNullSpaceDestroy(&pcbddc->NullSpace);CHKERRQ(ierr); 28 ierr = ISDestroy(&pcbddc->NeumannBoundaries);CHKERRQ(ierr); 29 ierr = ISDestroy(&pcbddc->DirichletBoundaries);CHKERRQ(ierr); 30 for (i=0;i<pcbddc->n_ISForDofs;i++) { 31 ierr = ISDestroy(&pcbddc->ISForDofs[i]);CHKERRQ(ierr); 32 } 33 ierr = PetscFree(pcbddc->ISForDofs);CHKERRQ(ierr); 34 PetscFunctionReturn(0); 35 } 36 37 #undef __FUNCT__ 38 #define __FUNCT__ "PCBDDCResetTopography" 39 PetscErrorCode PCBDDCResetTopography(PC pc) 40 { 41 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 42 PetscErrorCode ierr; 43 44 PetscFunctionBegin; 45 ierr = MatDestroy(&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 46 ierr = MatDestroy(&pcbddc->ConstraintMatrix);CHKERRQ(ierr); 47 ierr = PCBDDCGraphReset(pcbddc->mat_graph);CHKERRQ(ierr); 48 PetscFunctionReturn(0); 49 } 50 51 #undef __FUNCT__ 52 #define __FUNCT__ "PCBDDCResetSolvers" 53 PetscErrorCode PCBDDCResetSolvers(PC pc) 54 { 55 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 56 PetscErrorCode ierr; 57 58 PetscFunctionBegin; 59 ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr); 60 ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr); 61 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 62 ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr); 63 ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr); 64 ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr); 65 ierr = MatDestroy(&pcbddc->coarse_psi_B);CHKERRQ(ierr); 66 ierr = MatDestroy(&pcbddc->coarse_psi_D);CHKERRQ(ierr); 67 ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr); 68 ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr); 69 ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr); 70 ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr); 71 ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr); 72 ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr); 73 ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr); 74 ierr = ISDestroy(&pcbddc->is_R_local);CHKERRQ(ierr); 75 ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr); 76 ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr); 77 ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 78 ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr); 79 ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 80 ierr = PetscFree(pcbddc->replicated_local_primal_values);CHKERRQ(ierr); 81 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 82 ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr); 83 PetscFunctionReturn(0); 84 } 85 86 #undef __FUNCT__ 87 #define __FUNCT__ "PCBDDCCreateWorkVectors" 88 PetscErrorCode PCBDDCCreateWorkVectors(PC pc) 89 { 90 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 91 PC_IS *pcis = (PC_IS*)pc->data; 92 VecType impVecType; 93 PetscInt n_vertices,n_constraints,local_primal_size,n_R; 94 PetscErrorCode ierr; 95 96 PetscFunctionBegin; 97 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,NULL);CHKERRQ(ierr); 98 ierr = PCBDDCGetPrimalConstraintsLocalIdx(pc,&n_constraints,NULL);CHKERRQ(ierr); 99 local_primal_size = n_constraints+n_vertices; 100 n_R = pcis->n-n_vertices; 101 /* local work vectors */ 102 ierr = VecGetType(pcis->vec1_N,&impVecType);CHKERRQ(ierr); 103 ierr = VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);CHKERRQ(ierr); 104 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_R);CHKERRQ(ierr); 105 ierr = VecSetSizes(pcbddc->vec1_R,PETSC_DECIDE,n_R);CHKERRQ(ierr); 106 ierr = VecSetType(pcbddc->vec1_R,impVecType);CHKERRQ(ierr); 107 ierr = VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);CHKERRQ(ierr); 108 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_P);CHKERRQ(ierr); 109 ierr = VecSetSizes(pcbddc->vec1_P,PETSC_DECIDE,local_primal_size);CHKERRQ(ierr); 110 ierr = VecSetType(pcbddc->vec1_P,impVecType);CHKERRQ(ierr); 111 if (n_constraints) { 112 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_C);CHKERRQ(ierr); 113 ierr = VecSetSizes(pcbddc->vec1_C,PETSC_DECIDE,n_constraints);CHKERRQ(ierr); 114 ierr = VecSetType(pcbddc->vec1_C,impVecType);CHKERRQ(ierr); 115 } 116 PetscFunctionReturn(0); 117 } 118 119 #undef __FUNCT__ 120 #define __FUNCT__ "PCBDDCSetUpCoarseLocal" 121 PetscErrorCode PCBDDCSetUpCoarseLocal(PC pc) 122 { 123 PetscErrorCode ierr; 124 /* pointers to pcis and pcbddc */ 125 PC_IS* pcis = (PC_IS*)pc->data; 126 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 127 /* submatrices of local problem */ 128 Mat A_RV,A_VR,A_VV; 129 /* working matrices */ 130 Mat M1,M2,M3,C_CR; 131 /* working vectors */ 132 Vec vec1_C,vec2_C,vec1_V,vec2_V; 133 /* additional working stuff */ 134 IS is_aux; 135 ISLocalToGlobalMapping BtoNmap; 136 PetscScalar *coarse_submat_vals; /* TODO: use a PETSc matrix */ 137 const PetscScalar *array,*row_cmat_values; 138 const PetscInt *row_cmat_indices,*idx_R_local; 139 PetscInt *vertices,*idx_V_B,*auxindices; 140 PetscInt n_vertices,n_constraints,size_of_constraint; 141 PetscInt i,j,n_R,n_D,n_B; 142 PetscBool setsym=PETSC_FALSE,issym=PETSC_FALSE; 143 /* Vector and matrix types */ 144 VecType impVecType; 145 MatType impMatType; 146 /* some shortcuts to scalars */ 147 PetscScalar zero=0.0,one=1.0,m_one=-1.0; 148 /* for debugging purposes */ 149 PetscReal *coarsefunctions_errors,*constraints_errors; 150 151 PetscFunctionBegin; 152 /* get number of vertices and their local indices */ 153 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 154 n_constraints = pcbddc->local_primal_size-n_vertices; 155 /* Set Non-overlapping dimensions */ 156 n_B = pcis->n_B; n_D = pcis->n - n_B; 157 n_R = pcis->n-n_vertices; 158 159 /* Set types for local objects needed by BDDC precondtioner */ 160 impMatType = MATSEQDENSE; 161 ierr = VecGetType(pcis->vec1_N,&impVecType);CHKERRQ(ierr); 162 163 /* Allocating some extra storage just to be safe */ 164 ierr = PetscMalloc (pcis->n*sizeof(PetscInt),&auxindices);CHKERRQ(ierr); 165 for (i=0;i<pcis->n;i++) auxindices[i]=i; 166 167 /* vertices in boundary numbering */ 168 ierr = PetscMalloc(n_vertices*sizeof(PetscInt),&idx_V_B);CHKERRQ(ierr); 169 ierr = ISLocalToGlobalMappingCreateIS(pcis->is_B_local,&BtoNmap);CHKERRQ(ierr); 170 ierr = ISGlobalToLocalMappingApply(BtoNmap,IS_GTOLM_DROP,n_vertices,vertices,&i,idx_V_B);CHKERRQ(ierr); 171 ierr = ISLocalToGlobalMappingDestroy(&BtoNmap);CHKERRQ(ierr); 172 if (i != n_vertices) { 173 SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_SUP,"Error in boundary numbering for BDDC vertices! %d != %d\n",n_vertices,i); 174 } 175 176 /* some work vectors on vertices and/or constraints */ 177 if (n_vertices) { 178 ierr = VecCreate(PETSC_COMM_SELF,&vec1_V);CHKERRQ(ierr); 179 ierr = VecSetSizes(vec1_V,n_vertices,n_vertices);CHKERRQ(ierr); 180 ierr = VecSetType(vec1_V,impVecType);CHKERRQ(ierr); 181 ierr = VecDuplicate(vec1_V,&vec2_V);CHKERRQ(ierr); 182 } 183 if (n_constraints) { 184 ierr = VecDuplicate(pcbddc->vec1_C,&vec1_C);CHKERRQ(ierr); 185 ierr = VecDuplicate(pcbddc->vec1_C,&vec2_C);CHKERRQ(ierr); 186 } 187 188 /* Precompute stuffs needed for preprocessing and application of BDDC*/ 189 if (n_constraints) { 190 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->local_auxmat2);CHKERRQ(ierr); 191 ierr = MatSetSizes(pcbddc->local_auxmat2,n_R,n_constraints,PETSC_DECIDE,PETSC_DECIDE);CHKERRQ(ierr); 192 ierr = MatSetType(pcbddc->local_auxmat2,impMatType);CHKERRQ(ierr); 193 ierr = MatSetUp(pcbddc->local_auxmat2);CHKERRQ(ierr); 194 195 /* Extract constraints on R nodes: C_{CR} */ 196 ierr = ISCreateStride(PETSC_COMM_SELF,n_constraints,n_vertices,1,&is_aux);CHKERRQ(ierr); 197 ierr = MatGetSubMatrix(pcbddc->ConstraintMatrix,is_aux,pcbddc->is_R_local,MAT_INITIAL_MATRIX,&C_CR);CHKERRQ(ierr); 198 ierr = ISDestroy(&is_aux);CHKERRQ(ierr); 199 200 /* Assemble local_auxmat2 = - A_{RR}^{-1} C^T_{CR} needed by BDDC application */ 201 for (i=0;i<n_constraints;i++) { 202 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 203 /* Get row of constraint matrix in R numbering */ 204 ierr = MatGetRow(C_CR,i,&size_of_constraint,&row_cmat_indices,&row_cmat_values);CHKERRQ(ierr); 205 ierr = VecSetValues(pcbddc->vec1_R,size_of_constraint,row_cmat_indices,row_cmat_values,INSERT_VALUES);CHKERRQ(ierr); 206 ierr = MatRestoreRow(C_CR,i,&size_of_constraint,&row_cmat_indices,&row_cmat_values);CHKERRQ(ierr); 207 ierr = VecAssemblyBegin(pcbddc->vec1_R);CHKERRQ(ierr); 208 ierr = VecAssemblyEnd(pcbddc->vec1_R);CHKERRQ(ierr); 209 /* Solve for row of constraint matrix in R numbering */ 210 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 211 /* Set values in local_auxmat2 */ 212 ierr = VecGetArrayRead(pcbddc->vec2_R,&array);CHKERRQ(ierr); 213 ierr = MatSetValues(pcbddc->local_auxmat2,n_R,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 214 ierr = VecRestoreArrayRead(pcbddc->vec2_R,&array);CHKERRQ(ierr); 215 } 216 ierr = MatAssemblyBegin(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 217 ierr = MatAssemblyEnd(pcbddc->local_auxmat2,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 218 ierr = MatScale(pcbddc->local_auxmat2,m_one);CHKERRQ(ierr); 219 220 /* Assemble explicitly M1 = ( C_{CR} A_{RR}^{-1} C^T_{CR} )^{-1} needed in preproc */ 221 ierr = MatMatMult(C_CR,pcbddc->local_auxmat2,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&M3);CHKERRQ(ierr); 222 ierr = MatLUFactor(M3,NULL,NULL,NULL);CHKERRQ(ierr); 223 ierr = MatCreate(PETSC_COMM_SELF,&M1);CHKERRQ(ierr); 224 ierr = MatSetSizes(M1,n_constraints,n_constraints,n_constraints,n_constraints);CHKERRQ(ierr); 225 ierr = MatSetType(M1,impMatType);CHKERRQ(ierr); 226 ierr = MatSetUp(M1);CHKERRQ(ierr); 227 ierr = MatDuplicate(M1,MAT_DO_NOT_COPY_VALUES,&M2);CHKERRQ(ierr); 228 ierr = MatZeroEntries(M2);CHKERRQ(ierr); 229 ierr = VecSet(vec1_C,m_one);CHKERRQ(ierr); 230 ierr = MatDiagonalSet(M2,vec1_C,INSERT_VALUES);CHKERRQ(ierr); 231 ierr = MatMatSolve(M3,M2,M1);CHKERRQ(ierr); 232 ierr = MatDestroy(&M2);CHKERRQ(ierr); 233 ierr = MatDestroy(&M3);CHKERRQ(ierr); 234 /* Assemble local_auxmat1 = M1*C_{CR} needed by BDDC application in KSP and in preproc */ 235 ierr = MatMatMult(M1,C_CR,MAT_INITIAL_MATRIX,PETSC_DEFAULT,&pcbddc->local_auxmat1);CHKERRQ(ierr); 236 } 237 238 /* Get submatrices from subdomain matrix */ 239 if (n_vertices) { 240 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_vertices,vertices,PETSC_COPY_VALUES,&is_aux);CHKERRQ(ierr); 241 ierr = MatGetSubMatrix(pcbddc->local_mat,pcbddc->is_R_local,is_aux,MAT_INITIAL_MATRIX,&A_RV);CHKERRQ(ierr); 242 ierr = MatGetSubMatrix(pcbddc->local_mat,is_aux,pcbddc->is_R_local,MAT_INITIAL_MATRIX,&A_VR);CHKERRQ(ierr); 243 ierr = MatGetSubMatrix(pcbddc->local_mat,is_aux,is_aux,MAT_INITIAL_MATRIX,&A_VV);CHKERRQ(ierr); 244 ierr = ISDestroy(&is_aux);CHKERRQ(ierr); 245 } 246 247 /* Matrix of coarse basis functions (local) */ 248 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_B);CHKERRQ(ierr); 249 ierr = MatSetSizes(pcbddc->coarse_phi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 250 ierr = MatSetType(pcbddc->coarse_phi_B,impMatType);CHKERRQ(ierr); 251 ierr = MatSetUp(pcbddc->coarse_phi_B);CHKERRQ(ierr); 252 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag) { 253 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_phi_D);CHKERRQ(ierr); 254 ierr = MatSetSizes(pcbddc->coarse_phi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 255 ierr = MatSetType(pcbddc->coarse_phi_D,impMatType);CHKERRQ(ierr); 256 ierr = MatSetUp(pcbddc->coarse_phi_D);CHKERRQ(ierr); 257 } 258 259 if (pcbddc->dbg_flag) { 260 ierr = ISGetIndices(pcbddc->is_R_local,&idx_R_local);CHKERRQ(ierr); 261 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*coarsefunctions_errors),&coarsefunctions_errors);CHKERRQ(ierr); 262 ierr = PetscMalloc(2*pcbddc->local_primal_size*sizeof(*constraints_errors),&constraints_errors);CHKERRQ(ierr); 263 } 264 /* Subdomain contribution (Non-overlapping) to coarse matrix */ 265 ierr = PetscMalloc((pcbddc->local_primal_size)*(pcbddc->local_primal_size)*sizeof(PetscScalar),&coarse_submat_vals);CHKERRQ(ierr); 266 267 /* We are now ready to evaluate coarse basis functions and subdomain contribution to coarse problem */ 268 269 /* vertices */ 270 for (i=0;i<n_vertices;i++) { 271 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 272 ierr = VecSetValue(vec1_V,i,one,INSERT_VALUES);CHKERRQ(ierr); 273 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 274 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 275 /* simplified solution of saddle point problem with null rhs on constraints multipliers */ 276 ierr = MatMult(A_RV,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 277 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 278 ierr = VecScale(pcbddc->vec1_R,m_one);CHKERRQ(ierr); 279 if (n_constraints) { 280 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec1_R,vec1_C);CHKERRQ(ierr); 281 ierr = MatMultAdd(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 282 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 283 } 284 ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); 285 ierr = MatMultAdd(A_VV,vec1_V,vec2_V,vec2_V);CHKERRQ(ierr); 286 287 /* Set values in coarse basis function and subdomain part of coarse_mat */ 288 /* coarse basis functions */ 289 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 290 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 291 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 292 ierr = VecGetArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 293 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 294 ierr = VecRestoreArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 295 ierr = MatSetValue(pcbddc->coarse_phi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 296 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag) { 297 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 298 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 299 ierr = VecGetArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 300 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 301 ierr = VecRestoreArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 302 } 303 /* subdomain contribution to coarse matrix. WARNING -> column major ordering */ 304 ierr = VecGetArrayRead(vec2_V,&array);CHKERRQ(ierr); 305 ierr = PetscMemcpy(&coarse_submat_vals[i*pcbddc->local_primal_size],array,n_vertices*sizeof(PetscScalar));CHKERRQ(ierr); 306 ierr = VecRestoreArrayRead(vec2_V,&array);CHKERRQ(ierr); 307 if (n_constraints) { 308 ierr = VecGetArrayRead(vec1_C,&array);CHKERRQ(ierr); 309 ierr = PetscMemcpy(&coarse_submat_vals[i*pcbddc->local_primal_size+n_vertices],array,n_constraints*sizeof(PetscScalar));CHKERRQ(ierr); 310 ierr = VecRestoreArrayRead(vec1_C,&array);CHKERRQ(ierr); 311 } 312 313 /* check */ 314 if (pcbddc->dbg_flag) { 315 /* assemble subdomain vector on local nodes */ 316 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 317 ierr = VecGetArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 318 ierr = VecSetValues(pcis->vec1_N,n_R,idx_R_local,array,INSERT_VALUES);CHKERRQ(ierr); 319 ierr = VecRestoreArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 320 ierr = VecSetValue(pcis->vec1_N,vertices[i],one,INSERT_VALUES);CHKERRQ(ierr); 321 ierr = VecAssemblyBegin(pcis->vec1_N);CHKERRQ(ierr); 322 ierr = VecAssemblyEnd(pcis->vec1_N);CHKERRQ(ierr); 323 /* assemble subdomain vector of lagrange multipliers (i.e. primal nodes) */ 324 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 325 ierr = VecGetArrayRead(vec2_V,&array);CHKERRQ(ierr); 326 ierr = VecSetValues(pcbddc->vec1_P,n_vertices,auxindices,array,INSERT_VALUES);CHKERRQ(ierr); 327 ierr = VecRestoreArrayRead(vec2_V,&array);CHKERRQ(ierr); 328 if (n_constraints) { 329 ierr = VecGetArrayRead(vec1_C,&array);CHKERRQ(ierr); 330 ierr = VecSetValues(pcbddc->vec1_P,n_constraints,&auxindices[n_vertices],array,INSERT_VALUES);CHKERRQ(ierr); 331 ierr = VecRestoreArrayRead(vec1_C,&array);CHKERRQ(ierr); 332 } 333 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 334 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 335 ierr = VecScale(pcbddc->vec1_P,m_one);CHKERRQ(ierr); 336 /* check saddle point solution */ 337 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 338 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 339 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i]);CHKERRQ(ierr); 340 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 341 /* shift by the identity matrix */ 342 ierr = VecSetValue(pcbddc->vec1_P,i,m_one,ADD_VALUES);CHKERRQ(ierr); 343 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 344 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 345 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i]);CHKERRQ(ierr); 346 } 347 } 348 349 /* constraints */ 350 for (i=0;i<n_constraints;i++) { 351 ierr = VecSet(vec2_C,zero);CHKERRQ(ierr); 352 ierr = VecSetValue(vec2_C,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 353 ierr = VecAssemblyBegin(vec2_C);CHKERRQ(ierr); 354 ierr = VecAssemblyEnd(vec2_C);CHKERRQ(ierr); 355 /* simplified solution of saddle point problem with null rhs on vertices multipliers */ 356 ierr = MatMult(M1,vec2_C,vec1_C);CHKERRQ(ierr); 357 ierr = MatMult(pcbddc->local_auxmat2,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 358 ierr = VecScale(vec1_C,m_one);CHKERRQ(ierr); 359 if (n_vertices) { 360 ierr = MatMult(A_VR,pcbddc->vec1_R,vec2_V);CHKERRQ(ierr); 361 } 362 /* Set values in coarse basis function and subdomain part of coarse_mat */ 363 /* coarse basis functions */ 364 j = i+n_vertices; /* don't touch this! */ 365 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 366 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 367 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 368 ierr = VecGetArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 369 ierr = MatSetValues(pcbddc->coarse_phi_B,n_B,auxindices,1,&j,array,INSERT_VALUES);CHKERRQ(ierr); 370 ierr = VecRestoreArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 371 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag) { 372 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 373 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 374 ierr = VecGetArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 375 ierr = MatSetValues(pcbddc->coarse_phi_D,n_D,auxindices,1,&j,array,INSERT_VALUES);CHKERRQ(ierr); 376 ierr = VecRestoreArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 377 } 378 /* subdomain contribution to coarse matrix. WARNING -> column major ordering */ 379 if (n_vertices) { 380 ierr = VecGetArrayRead(vec2_V,&array);CHKERRQ(ierr); 381 ierr = PetscMemcpy(&coarse_submat_vals[j*pcbddc->local_primal_size],array,n_vertices*sizeof(PetscScalar));CHKERRQ(ierr); 382 ierr = VecRestoreArrayRead(vec2_V,&array);CHKERRQ(ierr); 383 } 384 ierr = VecGetArrayRead(vec1_C,&array);CHKERRQ(ierr); 385 ierr = PetscMemcpy(&coarse_submat_vals[j*pcbddc->local_primal_size+n_vertices],array,n_constraints*sizeof(PetscScalar));CHKERRQ(ierr); 386 ierr = VecRestoreArrayRead(vec1_C,&array);CHKERRQ(ierr); 387 388 if (pcbddc->dbg_flag) { 389 /* assemble subdomain vector on nodes */ 390 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 391 ierr = VecGetArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 392 ierr = VecSetValues(pcis->vec1_N,n_R,idx_R_local,array,INSERT_VALUES);CHKERRQ(ierr); 393 ierr = VecRestoreArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 394 ierr = VecAssemblyBegin(pcis->vec1_N);CHKERRQ(ierr); 395 ierr = VecAssemblyEnd(pcis->vec1_N);CHKERRQ(ierr); 396 /* assemble subdomain vector of lagrange multipliers */ 397 ierr = VecSet(pcbddc->vec1_P,zero);CHKERRQ(ierr); 398 if (n_vertices) { 399 ierr = VecGetArrayRead(vec2_V,&array);CHKERRQ(ierr); 400 ierr = VecSetValues(pcbddc->vec1_P,n_vertices,auxindices,array,INSERT_VALUES);CHKERRQ(ierr); 401 ierr = VecRestoreArrayRead(vec2_V,&array);CHKERRQ(ierr); 402 } 403 ierr = VecGetArrayRead(vec1_C,&array);CHKERRQ(ierr); 404 ierr = VecSetValues(pcbddc->vec1_P,n_constraints,&auxindices[n_vertices],array,INSERT_VALUES);CHKERRQ(ierr); 405 ierr = VecRestoreArrayRead(vec1_C,&array);CHKERRQ(ierr); 406 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 407 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 408 ierr = VecScale(pcbddc->vec1_P,m_one);CHKERRQ(ierr); 409 /* check saddle point solution */ 410 ierr = MatMult(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 411 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 412 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[j]);CHKERRQ(ierr); 413 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 414 /* shift by the identity matrix */ 415 ierr = VecSetValue(pcbddc->vec1_P,j,m_one,ADD_VALUES);CHKERRQ(ierr); 416 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 417 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 418 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[j]);CHKERRQ(ierr); 419 } 420 } 421 /* call assembling routines for local coarse basis */ 422 ierr = MatAssemblyBegin(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 423 ierr = MatAssemblyEnd(pcbddc->coarse_phi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 424 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag) { 425 ierr = MatAssemblyBegin(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 426 ierr = MatAssemblyEnd(pcbddc->coarse_phi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 427 } 428 429 /* compute other basis functions for non-symmetric problems */ 430 ierr = MatIsSymmetricKnown(pc->pmat,&setsym,&issym);CHKERRQ(ierr); 431 if (!setsym || (setsym && !issym)) { 432 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_B);CHKERRQ(ierr); 433 ierr = MatSetSizes(pcbddc->coarse_psi_B,n_B,pcbddc->local_primal_size,n_B,pcbddc->local_primal_size);CHKERRQ(ierr); 434 ierr = MatSetType(pcbddc->coarse_psi_B,impMatType);CHKERRQ(ierr); 435 ierr = MatSetUp(pcbddc->coarse_psi_B);CHKERRQ(ierr); 436 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag ) { 437 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->coarse_psi_D);CHKERRQ(ierr); 438 ierr = MatSetSizes(pcbddc->coarse_psi_D,n_D,pcbddc->local_primal_size,n_D,pcbddc->local_primal_size);CHKERRQ(ierr); 439 ierr = MatSetType(pcbddc->coarse_psi_D,impMatType);CHKERRQ(ierr); 440 ierr = MatSetUp(pcbddc->coarse_psi_D);CHKERRQ(ierr); 441 } 442 for (i=0;i<pcbddc->local_primal_size;i++) { 443 if (n_constraints) { 444 ierr = VecSet(vec1_C,zero);CHKERRQ(ierr); 445 for (j=0;j<n_constraints;j++) { 446 ierr = VecSetValue(vec1_C,j,coarse_submat_vals[(j+n_vertices)*pcbddc->local_primal_size+i],INSERT_VALUES);CHKERRQ(ierr); 447 } 448 ierr = VecAssemblyBegin(vec1_C);CHKERRQ(ierr); 449 ierr = VecAssemblyEnd(vec1_C);CHKERRQ(ierr); 450 } 451 if (i<n_vertices) { 452 ierr = VecSet(vec1_V,zero);CHKERRQ(ierr); 453 ierr = VecSetValue(vec1_V,i,m_one,INSERT_VALUES);CHKERRQ(ierr); 454 ierr = VecAssemblyBegin(vec1_V);CHKERRQ(ierr); 455 ierr = VecAssemblyEnd(vec1_V);CHKERRQ(ierr); 456 ierr = MatMultTranspose(A_VR,vec1_V,pcbddc->vec1_R);CHKERRQ(ierr); 457 if (n_constraints) { 458 ierr = MatMultTransposeAdd(C_CR,vec1_C,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 459 } 460 } else { 461 ierr = MatMultTranspose(C_CR,vec1_C,pcbddc->vec1_R);CHKERRQ(ierr); 462 } 463 ierr = KSPSolveTranspose(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec1_R);CHKERRQ(ierr); 464 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 465 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 466 ierr = VecScatterEnd(pcbddc->R_to_B,pcbddc->vec1_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 467 ierr = VecGetArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 468 ierr = MatSetValues(pcbddc->coarse_psi_B,n_B,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 469 ierr = VecRestoreArrayRead(pcis->vec1_B,&array);CHKERRQ(ierr); 470 if (i<n_vertices) { 471 ierr = MatSetValue(pcbddc->coarse_psi_B,idx_V_B[i],i,one,INSERT_VALUES);CHKERRQ(ierr); 472 } 473 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag) { 474 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 475 ierr = VecScatterEnd(pcbddc->R_to_D,pcbddc->vec1_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 476 ierr = VecGetArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 477 ierr = MatSetValues(pcbddc->coarse_psi_D,n_D,auxindices,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 478 ierr = VecRestoreArrayRead(pcis->vec1_D,&array);CHKERRQ(ierr); 479 } 480 481 if (pcbddc->dbg_flag) { 482 /* assemble subdomain vector on nodes */ 483 ierr = VecSet(pcis->vec1_N,zero);CHKERRQ(ierr); 484 ierr = VecGetArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 485 ierr = VecSetValues(pcis->vec1_N,n_R,idx_R_local,array,INSERT_VALUES);CHKERRQ(ierr); 486 ierr = VecRestoreArrayRead(pcbddc->vec1_R,&array);CHKERRQ(ierr); 487 if (i<n_vertices) { 488 ierr = VecSetValue(pcis->vec1_N,vertices[i],one,INSERT_VALUES);CHKERRQ(ierr); 489 } 490 ierr = VecAssemblyBegin(pcis->vec1_N);CHKERRQ(ierr); 491 ierr = VecAssemblyEnd(pcis->vec1_N);CHKERRQ(ierr); 492 /* assemble subdomain vector of lagrange multipliers */ 493 for (j=0;j<pcbddc->local_primal_size;j++) { 494 ierr = VecSetValue(pcbddc->vec1_P,j,-coarse_submat_vals[j*pcbddc->local_primal_size+i],INSERT_VALUES);CHKERRQ(ierr); 495 } 496 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 497 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 498 /* check saddle point solution */ 499 ierr = MatMultTranspose(pcbddc->local_mat,pcis->vec1_N,pcis->vec2_N);CHKERRQ(ierr); 500 ierr = MatMultTransposeAdd(pcbddc->ConstraintMatrix,pcbddc->vec1_P,pcis->vec2_N,pcis->vec2_N);CHKERRQ(ierr); 501 ierr = VecNorm(pcis->vec2_N,NORM_INFINITY,&coarsefunctions_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 502 ierr = MatMult(pcbddc->ConstraintMatrix,pcis->vec1_N,pcbddc->vec1_P);CHKERRQ(ierr); 503 /* shift by the identity matrix */ 504 ierr = VecSetValue(pcbddc->vec1_P,i,m_one,ADD_VALUES);CHKERRQ(ierr); 505 ierr = VecAssemblyBegin(pcbddc->vec1_P);CHKERRQ(ierr); 506 ierr = VecAssemblyEnd(pcbddc->vec1_P);CHKERRQ(ierr); 507 ierr = VecNorm(pcbddc->vec1_P,NORM_INFINITY,&constraints_errors[i+pcbddc->local_primal_size]);CHKERRQ(ierr); 508 } 509 } 510 ierr = MatAssemblyBegin(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 511 ierr = MatAssemblyEnd(pcbddc->coarse_psi_B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 512 if ( pcbddc->inexact_prec_type || pcbddc->dbg_flag ) { 513 ierr = MatAssemblyBegin(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 514 ierr = MatAssemblyEnd(pcbddc->coarse_psi_D,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 515 } 516 } 517 ierr = PetscFree(idx_V_B);CHKERRQ(ierr); 518 /* Checking coarse_sub_mat and coarse basis functios */ 519 /* Symmetric case : It should be \Phi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 520 /* Non-symmetric case : It should be \Psi^{(j)^T} A^{(j)} \Phi^{(j)}=coarse_sub_mat */ 521 if (pcbddc->dbg_flag) { 522 Mat coarse_sub_mat; 523 Mat AUXMAT,TM1,TM2,TM3,TM4; 524 Mat coarse_phi_D,coarse_phi_B; 525 Mat coarse_psi_D,coarse_psi_B; 526 Mat A_II,A_BB,A_IB,A_BI; 527 MatType checkmattype=MATSEQAIJ; 528 PetscReal real_value; 529 530 ierr = MatConvert(pcis->A_II,checkmattype,MAT_INITIAL_MATRIX,&A_II);CHKERRQ(ierr); 531 ierr = MatConvert(pcis->A_IB,checkmattype,MAT_INITIAL_MATRIX,&A_IB);CHKERRQ(ierr); 532 ierr = MatConvert(pcis->A_BI,checkmattype,MAT_INITIAL_MATRIX,&A_BI);CHKERRQ(ierr); 533 ierr = MatConvert(pcis->A_BB,checkmattype,MAT_INITIAL_MATRIX,&A_BB);CHKERRQ(ierr); 534 ierr = MatConvert(pcbddc->coarse_phi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_D);CHKERRQ(ierr); 535 ierr = MatConvert(pcbddc->coarse_phi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_phi_B);CHKERRQ(ierr); 536 if (pcbddc->coarse_psi_B) { 537 ierr = MatConvert(pcbddc->coarse_psi_D,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_D);CHKERRQ(ierr); 538 ierr = MatConvert(pcbddc->coarse_psi_B,checkmattype,MAT_INITIAL_MATRIX,&coarse_psi_B);CHKERRQ(ierr); 539 } 540 ierr = MatCreateSeqDense(PETSC_COMM_SELF,pcbddc->local_primal_size,pcbddc->local_primal_size,coarse_submat_vals,&coarse_sub_mat);CHKERRQ(ierr); 541 542 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 543 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Check coarse sub mat and local basis functions\n");CHKERRQ(ierr); 544 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 545 if (pcbddc->coarse_psi_B) { 546 ierr = MatMatMult(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 547 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 548 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 549 ierr = MatMatMult(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 550 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 551 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 552 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 553 ierr = MatTransposeMatMult(coarse_psi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 554 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 555 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 556 ierr = MatTransposeMatMult(coarse_psi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 557 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 558 } else { 559 ierr = MatPtAP(A_II,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&TM1);CHKERRQ(ierr); 560 ierr = MatPtAP(A_BB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&TM2);CHKERRQ(ierr); 561 ierr = MatMatMult(A_IB,coarse_phi_B,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 562 ierr = MatTransposeMatMult(coarse_phi_D,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM3);CHKERRQ(ierr); 563 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 564 ierr = MatMatMult(A_BI,coarse_phi_D,MAT_INITIAL_MATRIX,1.0,&AUXMAT);CHKERRQ(ierr); 565 ierr = MatTransposeMatMult(coarse_phi_B,AUXMAT,MAT_INITIAL_MATRIX,1.0,&TM4);CHKERRQ(ierr); 566 ierr = MatDestroy(&AUXMAT);CHKERRQ(ierr); 567 } 568 ierr = MatAXPY(TM1,one,TM2,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 569 ierr = MatAXPY(TM1,one,TM3,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 570 ierr = MatAXPY(TM1,one,TM4,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr); 571 ierr = MatConvert(TM1,MATSEQDENSE,MAT_REUSE_MATRIX,&TM1);CHKERRQ(ierr); 572 ierr = MatAXPY(TM1,m_one,coarse_sub_mat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 573 ierr = MatNorm(TM1,NORM_INFINITY,&real_value);CHKERRQ(ierr); 574 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"----------------------------------\n");CHKERRQ(ierr); 575 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d \n",PetscGlobalRank);CHKERRQ(ierr); 576 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"matrix error = % 1.14e\n",real_value);CHKERRQ(ierr); 577 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"coarse functions (phi) errors\n");CHKERRQ(ierr); 578 for (i=0;i<pcbddc->local_primal_size;i++) { 579 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local %02d-th function error = % 1.14e\n",i,coarsefunctions_errors[i]);CHKERRQ(ierr); 580 } 581 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"constraints (phi) errors\n");CHKERRQ(ierr); 582 for (i=0;i<pcbddc->local_primal_size;i++) { 583 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local %02d-th function error = % 1.14e\n",i,constraints_errors[i]);CHKERRQ(ierr); 584 } 585 if (pcbddc->coarse_psi_B) { 586 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"coarse functions (psi) errors\n");CHKERRQ(ierr); 587 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 588 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,coarsefunctions_errors[i]);CHKERRQ(ierr); 589 } 590 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"constraints (psi) errors\n");CHKERRQ(ierr); 591 for (i=pcbddc->local_primal_size;i<2*pcbddc->local_primal_size;i++) { 592 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local %02d-th function error = % 1.14e\n",i-pcbddc->local_primal_size,constraints_errors[i]);CHKERRQ(ierr); 593 } 594 } 595 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 596 ierr = MatDestroy(&A_II);CHKERRQ(ierr); 597 ierr = MatDestroy(&A_BB);CHKERRQ(ierr); 598 ierr = MatDestroy(&A_IB);CHKERRQ(ierr); 599 ierr = MatDestroy(&A_BI);CHKERRQ(ierr); 600 ierr = MatDestroy(&TM1);CHKERRQ(ierr); 601 ierr = MatDestroy(&TM2);CHKERRQ(ierr); 602 ierr = MatDestroy(&TM3);CHKERRQ(ierr); 603 ierr = MatDestroy(&TM4);CHKERRQ(ierr); 604 ierr = MatDestroy(&coarse_phi_D);CHKERRQ(ierr); 605 ierr = MatDestroy(&coarse_phi_B);CHKERRQ(ierr); 606 if (pcbddc->coarse_psi_B) { 607 ierr = MatDestroy(&coarse_psi_D);CHKERRQ(ierr); 608 ierr = MatDestroy(&coarse_psi_B);CHKERRQ(ierr); 609 } 610 ierr = MatDestroy(&coarse_sub_mat);CHKERRQ(ierr); 611 ierr = ISRestoreIndices(pcbddc->is_R_local,&idx_R_local);CHKERRQ(ierr); 612 ierr = PetscFree(coarsefunctions_errors);CHKERRQ(ierr); 613 ierr = PetscFree(constraints_errors);CHKERRQ(ierr); 614 } 615 /* free memory */ 616 if (n_vertices) { 617 ierr = PetscFree(vertices);CHKERRQ(ierr); 618 ierr = VecDestroy(&vec1_V);CHKERRQ(ierr); 619 ierr = VecDestroy(&vec2_V);CHKERRQ(ierr); 620 ierr = MatDestroy(&A_RV);CHKERRQ(ierr); 621 ierr = MatDestroy(&A_VR);CHKERRQ(ierr); 622 ierr = MatDestroy(&A_VV);CHKERRQ(ierr); 623 } 624 if (n_constraints) { 625 ierr = VecDestroy(&vec1_C);CHKERRQ(ierr); 626 ierr = VecDestroy(&vec2_C);CHKERRQ(ierr); 627 ierr = MatDestroy(&M1);CHKERRQ(ierr); 628 ierr = MatDestroy(&C_CR);CHKERRQ(ierr); 629 } 630 ierr = PetscFree(auxindices);CHKERRQ(ierr); 631 /* create coarse matrix and data structures for message passing associated actual choice of coarse problem type */ 632 ierr = PCBDDCSetUpCoarseEnvironment(pc,coarse_submat_vals);CHKERRQ(ierr); 633 ierr = PetscFree(coarse_submat_vals);CHKERRQ(ierr); 634 PetscFunctionReturn(0); 635 } 636 637 #undef __FUNCT__ 638 #define __FUNCT__ "PCBDDCSetUpLocalMatrices" 639 PetscErrorCode PCBDDCSetUpLocalMatrices(PC pc) 640 { 641 PC_IS* pcis = (PC_IS*)(pc->data); 642 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 643 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 644 /* manage repeated solves */ 645 MatReuse reuse; 646 MatStructure matstruct; 647 PetscErrorCode ierr; 648 649 PetscFunctionBegin; 650 /* get mat flags */ 651 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 652 reuse = MAT_INITIAL_MATRIX; 653 if (pc->setupcalled) { 654 /* when matstruct is SAME_PRECONDITIONER, we shouldn't be here */ 655 if (matstruct == SAME_PRECONDITIONER) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"This should not happen"); 656 if (matstruct == SAME_NONZERO_PATTERN) { 657 reuse = MAT_REUSE_MATRIX; 658 } else { 659 reuse = MAT_INITIAL_MATRIX; 660 } 661 } 662 if (reuse == MAT_INITIAL_MATRIX) { 663 ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr); 664 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 665 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 666 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 667 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 668 } 669 670 /* transform local matrices if needed */ 671 if (pcbddc->use_change_of_basis) { 672 Mat change_mat_all; 673 PetscScalar *row_cmat_values; 674 PetscInt *row_cmat_indices; 675 PetscInt *nnz,*is_indices,*temp_indices; 676 PetscInt i,j,k,n_D,n_B; 677 678 /* Get Non-overlapping dimensions */ 679 n_B = pcis->n_B; 680 n_D = pcis->n-n_B; 681 682 /* compute nonzero structure of change of basis on all local nodes */ 683 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 684 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 685 for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1; 686 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 687 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 688 k=1; 689 for (i=0;i<n_B;i++) { 690 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 691 nnz[is_indices[i]]=j; 692 if (k < j) k = j; 693 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 694 } 695 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 696 /* assemble change of basis matrix on the whole set of local dofs */ 697 ierr = PetscMalloc(k*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 698 ierr = MatCreate(PETSC_COMM_SELF,&change_mat_all);CHKERRQ(ierr); 699 ierr = MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);CHKERRQ(ierr); 700 ierr = MatSetType(change_mat_all,MATSEQAIJ);CHKERRQ(ierr); 701 ierr = MatSeqAIJSetPreallocation(change_mat_all,0,nnz);CHKERRQ(ierr); 702 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 703 for (i=0;i<n_D;i++) { 704 ierr = MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 705 } 706 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 707 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 708 for (i=0;i<n_B;i++) { 709 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 710 for (k=0; k<j; k++) temp_indices[k]=is_indices[row_cmat_indices[k]]; 711 ierr = MatSetValues(change_mat_all,1,&is_indices[i],j,temp_indices,row_cmat_values,INSERT_VALUES);CHKERRQ(ierr); 712 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 713 } 714 ierr = MatAssemblyBegin(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 715 ierr = MatAssemblyEnd(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 716 /* TODO: HOW TO WORK WITH BAIJ? PtAP not provided */ 717 ierr = MatGetBlockSize(matis->A,&i);CHKERRQ(ierr); 718 if (i==1) { 719 ierr = MatPtAP(matis->A,change_mat_all,reuse,2.0,&pcbddc->local_mat);CHKERRQ(ierr); 720 } else { 721 Mat work_mat; 722 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&work_mat);CHKERRQ(ierr); 723 ierr = MatPtAP(work_mat,change_mat_all,reuse,2.0,&pcbddc->local_mat);CHKERRQ(ierr); 724 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 725 } 726 ierr = MatDestroy(&change_mat_all);CHKERRQ(ierr); 727 ierr = PetscFree(nnz);CHKERRQ(ierr); 728 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 729 } else { 730 /* without change of basis, the local matrix is unchanged */ 731 if (!pcbddc->local_mat) { 732 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 733 pcbddc->local_mat = matis->A; 734 } 735 } 736 737 /* get submatrices */ 738 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);CHKERRQ(ierr); 739 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);CHKERRQ(ierr); 740 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);CHKERRQ(ierr); 741 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);CHKERRQ(ierr); 742 PetscFunctionReturn(0); 743 } 744 745 #undef __FUNCT__ 746 #define __FUNCT__ "PCBDDCSetUpLocalScatters" 747 PetscErrorCode PCBDDCSetUpLocalScatters(PC pc) 748 { 749 PC_IS* pcis = (PC_IS*)(pc->data); 750 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 751 IS is_aux1,is_aux2; 752 PetscInt *vertices,*aux_array1,*aux_array2,*is_indices,*idx_R_local; 753 PetscInt n_vertices,n_constraints,i,j,n_R,n_D,n_B; 754 PetscBool *array_bool; 755 PetscErrorCode ierr; 756 757 PetscFunctionBegin; 758 /* Set Non-overlapping dimensions */ 759 n_B = pcis->n_B; n_D = pcis->n - n_B; 760 /* get vertex indices from constraint matrix */ 761 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 762 /* Set number of constraints */ 763 n_constraints = pcbddc->local_primal_size-n_vertices; 764 /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */ 765 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&array_bool);CHKERRQ(ierr); 766 for (i=0;i<pcis->n;i++) array_bool[i] = PETSC_TRUE; 767 for (i=0;i<n_vertices;i++) array_bool[vertices[i]] = PETSC_FALSE; 768 ierr = PetscMalloc((pcis->n-n_vertices)*sizeof(PetscInt),&idx_R_local);CHKERRQ(ierr); 769 for (i=0, n_R=0; i<pcis->n; i++) { 770 if (array_bool[i]) { 771 idx_R_local[n_R] = i; 772 n_R++; 773 } 774 } 775 ierr = PetscFree(vertices);CHKERRQ(ierr); 776 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_OWN_POINTER,&pcbddc->is_R_local);CHKERRQ(ierr); 777 778 /* print some info if requested */ 779 if (pcbddc->dbg_flag) { 780 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 781 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 782 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);CHKERRQ(ierr); 783 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);CHKERRQ(ierr); 784 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"r_size = %d, v_size = %d, constraints = %d, local_primal_size = %d\n",n_R,n_vertices,n_constraints,pcbddc->local_primal_size);CHKERRQ(ierr); 785 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"pcbddc->n_vertices = %d, pcbddc->n_constraints = %d\n",pcbddc->n_vertices,pcbddc->n_constraints);CHKERRQ(ierr); 786 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 787 } 788 789 /* VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ 790 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 791 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array2);CHKERRQ(ierr); 792 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 793 for (i=0; i<n_D; i++) array_bool[is_indices[i]] = PETSC_FALSE; 794 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 795 for (i=0, j=0; i<n_R; i++) { 796 if (array_bool[idx_R_local[i]]) { 797 aux_array1[j] = i; 798 j++; 799 } 800 } 801 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_OWN_POINTER,&is_aux1);CHKERRQ(ierr); 802 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 803 for (i=0, j=0; i<n_B; i++) { 804 if (array_bool[is_indices[i]]) { 805 aux_array2[j] = i; j++; 806 } 807 } 808 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 809 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array2,PETSC_OWN_POINTER,&is_aux2);CHKERRQ(ierr); 810 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);CHKERRQ(ierr); 811 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 812 ierr = ISDestroy(&is_aux2);CHKERRQ(ierr); 813 814 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag ) { 815 ierr = PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 816 for (i=0, j=0; i<n_R; i++) { 817 if (!array_bool[idx_R_local[i]]) { 818 aux_array1[j] = i; 819 j++; 820 } 821 } 822 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_OWN_POINTER,&is_aux1);CHKERRQ(ierr); 823 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);CHKERRQ(ierr); 824 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 825 } 826 ierr = PetscFree(array_bool);CHKERRQ(ierr); 827 PetscFunctionReturn(0); 828 } 829 830 #undef __FUNCT__ 831 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 832 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use) 833 { 834 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 835 836 PetscFunctionBegin; 837 pcbddc->use_exact_dirichlet=use; 838 PetscFunctionReturn(0); 839 } 840 841 #undef __FUNCT__ 842 #define __FUNCT__ "PCBDDCSetUpLocalSolvers" 843 PetscErrorCode PCBDDCSetUpLocalSolvers(PC pc) 844 { 845 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 846 PC_IS *pcis = (PC_IS*)pc->data; 847 PC pc_temp; 848 Mat A_RR; 849 Vec vec1,vec2,vec3; 850 MatStructure matstruct; 851 PetscScalar m_one = -1.0; 852 PetscReal value; 853 PetscInt n_D,n_R,use_exact,use_exact_reduced; 854 PetscErrorCode ierr; 855 856 PetscFunctionBegin; 857 /* Creating PC contexts for local Dirichlet and Neumann problems */ 858 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 859 860 /* DIRICHLET PROBLEM */ 861 /* Matrix for Dirichlet problem is pcis->A_II */ 862 ierr = ISGetSize(pcis->is_I_local,&n_D);CHKERRQ(ierr); 863 if (!pcbddc->ksp_D) { /* create object if not yet build */ 864 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 865 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 866 /* default */ 867 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 868 ierr = KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");CHKERRQ(ierr); 869 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 870 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 871 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 872 } 873 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,matstruct);CHKERRQ(ierr); 874 /* Allow user's customization */ 875 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 876 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 877 if (!n_D) { 878 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 879 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 880 } 881 /* Set Up KSP for Dirichlet problem of BDDC */ 882 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 883 /* set ksp_D into pcis data */ 884 ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr); 885 ierr = PetscObjectReference((PetscObject)pcbddc->ksp_D);CHKERRQ(ierr); 886 pcis->ksp_D = pcbddc->ksp_D; 887 888 /* NEUMANN PROBLEM */ 889 /* Matrix for Neumann problem is A_RR -> we need to create it */ 890 ierr = ISGetSize(pcbddc->is_R_local,&n_R);CHKERRQ(ierr); 891 ierr = MatGetSubMatrix(pcbddc->local_mat,pcbddc->is_R_local,pcbddc->is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 892 if (!pcbddc->ksp_R) { /* create object if not yet build */ 893 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 894 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 895 /* default */ 896 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 897 ierr = KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");CHKERRQ(ierr); 898 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 899 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 900 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 901 } 902 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,matstruct);CHKERRQ(ierr); 903 /* Allow user's customization */ 904 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 905 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 906 if (!n_R) { 907 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 908 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 909 } 910 /* Set Up KSP for Neumann problem of BDDC */ 911 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 912 913 /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */ 914 915 /* Dirichlet */ 916 ierr = MatGetVecs(pcis->A_II,&vec1,&vec2);CHKERRQ(ierr); 917 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 918 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 919 ierr = MatMult(pcis->A_II,vec1,vec2);CHKERRQ(ierr); 920 ierr = KSPSolve(pcbddc->ksp_D,vec2,vec3);CHKERRQ(ierr); 921 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 922 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 923 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 924 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 925 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 926 /* need to be adapted? */ 927 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 928 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 929 ierr = PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);CHKERRQ(ierr); 930 /* print info */ 931 if (pcbddc->dbg_flag) { 932 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 933 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 934 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Checking solution of Dirichlet and Neumann problems\n");CHKERRQ(ierr); 935 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 936 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 937 } 938 if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) { 939 ierr = PCBDDCNullSpaceAssembleCorrection(pc,pcis->is_I_local);CHKERRQ(ierr); 940 } 941 942 /* Neumann */ 943 ierr = MatGetVecs(A_RR,&vec1,&vec2);CHKERRQ(ierr); 944 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 945 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 946 ierr = MatMult(A_RR,vec1,vec2);CHKERRQ(ierr); 947 ierr = KSPSolve(pcbddc->ksp_R,vec2,vec3);CHKERRQ(ierr); 948 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 949 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 950 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 951 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 952 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 953 /* need to be adapted? */ 954 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 955 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 956 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 957 /* print info */ 958 if (pcbddc->dbg_flag) { 959 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 960 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 961 } 962 if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */ 963 ierr = PCBDDCNullSpaceAssembleCorrection(pc,pcbddc->is_R_local);CHKERRQ(ierr); 964 } 965 966 /* free Neumann problem's matrix */ 967 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 968 PetscFunctionReturn(0); 969 } 970 971 #undef __FUNCT__ 972 #define __FUNCT__ "PCBDDCSolveSaddlePoint" 973 static PetscErrorCode PCBDDCSolveSaddlePoint(PC pc) 974 { 975 PetscErrorCode ierr; 976 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 977 978 PetscFunctionBegin; 979 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 980 if (pcbddc->local_auxmat1) { 981 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec2_R,pcbddc->vec1_C);CHKERRQ(ierr); 982 ierr = MatMultAdd(pcbddc->local_auxmat2,pcbddc->vec1_C,pcbddc->vec2_R,pcbddc->vec2_R);CHKERRQ(ierr); 983 } 984 PetscFunctionReturn(0); 985 } 986 987 #undef __FUNCT__ 988 #define __FUNCT__ "PCBDDCApplyInterfacePreconditioner" 989 PetscErrorCode PCBDDCApplyInterfacePreconditioner(PC pc) 990 { 991 PetscErrorCode ierr; 992 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 993 PC_IS* pcis = (PC_IS*) (pc->data); 994 const PetscScalar zero = 0.0; 995 996 PetscFunctionBegin; 997 /* Application of PHI^T (or PSI^T) */ 998 if (pcbddc->coarse_psi_B) { 999 ierr = MatMultTranspose(pcbddc->coarse_psi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 1000 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_psi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 1001 } else { 1002 ierr = MatMultTranspose(pcbddc->coarse_phi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 1003 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_phi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 1004 } 1005 /* Scatter data of coarse_rhs */ 1006 if (pcbddc->coarse_rhs) { ierr = VecSet(pcbddc->coarse_rhs,zero);CHKERRQ(ierr); } 1007 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1008 1009 /* Local solution on R nodes */ 1010 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 1011 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1012 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1013 if (pcbddc->inexact_prec_type) { 1014 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1015 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1016 } 1017 ierr = PCBDDCSolveSaddlePoint(pc);CHKERRQ(ierr); 1018 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 1019 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1020 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1021 if (pcbddc->inexact_prec_type) { 1022 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1023 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1024 } 1025 1026 /* Coarse solution */ 1027 ierr = PCBDDCScatterCoarseDataEnd(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1028 if (pcbddc->coarse_rhs) { /* TODO remove null space when doing multilevel */ 1029 ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 1030 } 1031 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1032 ierr = PCBDDCScatterCoarseDataEnd (pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1033 1034 /* Sum contributions from two levels */ 1035 ierr = MatMultAdd(pcbddc->coarse_phi_B,pcbddc->vec1_P,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 1036 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcbddc->coarse_phi_D,pcbddc->vec1_P,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 1037 PetscFunctionReturn(0); 1038 } 1039 1040 #undef __FUNCT__ 1041 #define __FUNCT__ "PCBDDCScatterCoarseDataBegin" 1042 PetscErrorCode PCBDDCScatterCoarseDataBegin(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 1043 { 1044 PetscErrorCode ierr; 1045 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 1046 1047 PetscFunctionBegin; 1048 switch (pcbddc->coarse_communications_type) { 1049 case SCATTERS_BDDC: 1050 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 1051 break; 1052 case GATHERS_BDDC: 1053 break; 1054 } 1055 PetscFunctionReturn(0); 1056 } 1057 1058 #undef __FUNCT__ 1059 #define __FUNCT__ "PCBDDCScatterCoarseDataEnd" 1060 PetscErrorCode PCBDDCScatterCoarseDataEnd(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 1061 { 1062 PetscErrorCode ierr; 1063 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 1064 PetscScalar* array_to; 1065 PetscScalar* array_from; 1066 MPI_Comm comm; 1067 PetscInt i; 1068 1069 PetscFunctionBegin; 1070 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 1071 switch (pcbddc->coarse_communications_type) { 1072 case SCATTERS_BDDC: 1073 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 1074 break; 1075 case GATHERS_BDDC: 1076 if (vec_from) { 1077 ierr = VecGetArray(vec_from,&array_from);CHKERRQ(ierr); 1078 } 1079 if (vec_to) { 1080 ierr = VecGetArray(vec_to,&array_to);CHKERRQ(ierr); 1081 } 1082 switch(pcbddc->coarse_problem_type){ 1083 case SEQUENTIAL_BDDC: 1084 if (smode == SCATTER_FORWARD) { 1085 ierr = MPI_Gatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 1086 if (vec_to) { 1087 if (imode == ADD_VALUES) { 1088 for (i=0;i<pcbddc->replicated_primal_size;i++) { 1089 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 1090 } 1091 } else { 1092 for (i=0;i<pcbddc->replicated_primal_size;i++) { 1093 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 1094 } 1095 } 1096 } 1097 } else { 1098 if (vec_from) { 1099 if (imode == ADD_VALUES) { 1100 MPI_Comm vec_from_comm; 1101 ierr = PetscObjectGetComm((PetscObject)(vec_from),&vec_from_comm);CHKERRQ(ierr); 1102 SETERRQ2(vec_from_comm,PETSC_ERR_SUP,"Unsupported insert mode ADD_VALUES for SCATTER_REVERSE in %s for case %d\n",__FUNCT__,pcbddc->coarse_problem_type); 1103 } 1104 for (i=0;i<pcbddc->replicated_primal_size;i++) { 1105 pcbddc->replicated_local_primal_values[i]=array_from[pcbddc->replicated_local_primal_indices[i]]; 1106 } 1107 } 1108 ierr = MPI_Scatterv(&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,&array_to[0],pcbddc->local_primal_size,MPIU_SCALAR,0,comm);CHKERRQ(ierr); 1109 } 1110 break; 1111 case REPLICATED_BDDC: 1112 if (smode == SCATTER_FORWARD) { 1113 ierr = MPI_Allgatherv(&array_from[0],pcbddc->local_primal_size,MPIU_SCALAR,&pcbddc->replicated_local_primal_values[0],pcbddc->local_primal_sizes,pcbddc->local_primal_displacements,MPIU_SCALAR,comm);CHKERRQ(ierr); 1114 if (imode == ADD_VALUES) { 1115 for (i=0;i<pcbddc->replicated_primal_size;i++) { 1116 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 1117 } 1118 } else { 1119 for (i=0;i<pcbddc->replicated_primal_size;i++) { 1120 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 1121 } 1122 } 1123 } else { /* no communications needed for SCATTER_REVERSE since needed data is already present */ 1124 if (imode == ADD_VALUES) { 1125 for (i=0;i<pcbddc->local_primal_size;i++) { 1126 array_to[i]+=array_from[pcbddc->local_primal_indices[i]]; 1127 } 1128 } else { 1129 for (i=0;i<pcbddc->local_primal_size;i++) { 1130 array_to[i]=array_from[pcbddc->local_primal_indices[i]]; 1131 } 1132 } 1133 } 1134 break; 1135 case MULTILEVEL_BDDC: 1136 break; 1137 case PARALLEL_BDDC: 1138 break; 1139 } 1140 if (vec_from) { 1141 ierr = VecRestoreArray(vec_from,&array_from);CHKERRQ(ierr); 1142 } 1143 if (vec_to) { 1144 ierr = VecRestoreArray(vec_to,&array_to);CHKERRQ(ierr); 1145 } 1146 break; 1147 } 1148 PetscFunctionReturn(0); 1149 } 1150 1151 /* uncomment for testing purposes */ 1152 /* #define PETSC_MISSING_LAPACK_GESVD 1 */ 1153 #undef __FUNCT__ 1154 #define __FUNCT__ "PCBDDCConstraintsSetUp" 1155 PetscErrorCode PCBDDCConstraintsSetUp(PC pc) 1156 { 1157 PetscErrorCode ierr; 1158 PC_IS* pcis = (PC_IS*)(pc->data); 1159 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 1160 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 1161 /* constraint and (optionally) change of basis matrix implemented as SeqAIJ */ 1162 MatType impMatType=MATSEQAIJ; 1163 /* one and zero */ 1164 PetscScalar one=1.0,zero=0.0; 1165 /* space to store constraints and their local indices */ 1166 PetscScalar *temp_quadrature_constraint; 1167 PetscInt *temp_indices,*temp_indices_to_constraint,*temp_indices_to_constraint_B; 1168 /* iterators */ 1169 PetscInt i,j,k,total_counts,temp_start_ptr; 1170 /* stuff to store connected components stored in pcbddc->mat_graph */ 1171 IS ISForVertices,*ISForFaces,*ISForEdges,*used_IS; 1172 PetscInt n_ISForFaces,n_ISForEdges; 1173 PetscBool get_faces,get_edges,get_vertices; 1174 /* near null space stuff */ 1175 MatNullSpace nearnullsp; 1176 const Vec *nearnullvecs; 1177 Vec *localnearnullsp; 1178 PetscBool nnsp_has_cnst; 1179 PetscInt nnsp_size; 1180 PetscScalar *array; 1181 /* BLAS integers */ 1182 PetscBLASInt lwork,lierr; 1183 PetscBLASInt Blas_N,Blas_M,Blas_K,Blas_one=1; 1184 PetscBLASInt Blas_LDA,Blas_LDB,Blas_LDC; 1185 /* LAPACK working arrays for SVD or POD */ 1186 PetscBool skip_lapack; 1187 PetscScalar *work; 1188 PetscReal *singular_vals; 1189 #if defined(PETSC_USE_COMPLEX) 1190 PetscReal *rwork; 1191 #endif 1192 #if defined(PETSC_MISSING_LAPACK_GESVD) 1193 PetscBLASInt Blas_one_2=1; 1194 PetscScalar *temp_basis,*correlation_mat; 1195 #else 1196 PetscBLASInt dummy_int_1=1,dummy_int_2=1; 1197 PetscScalar dummy_scalar_1=0.0,dummy_scalar_2=0.0; 1198 #endif 1199 /* change of basis */ 1200 PetscInt *aux_primal_numbering,*aux_primal_minloc,*global_indices; 1201 PetscBool boolforchange,*change_basis,*touched; 1202 /* auxiliary stuff */ 1203 PetscInt *nnz,*is_indices,*local_to_B; 1204 /* some quantities */ 1205 PetscInt n_vertices,total_primal_vertices; 1206 PetscInt size_of_constraint,max_size_of_constraint,max_constraints,temp_constraints; 1207 1208 1209 PetscFunctionBegin; 1210 /* Get index sets for faces, edges and vertices from graph */ 1211 get_faces = PETSC_TRUE; 1212 get_edges = PETSC_TRUE; 1213 get_vertices = PETSC_TRUE; 1214 if (pcbddc->vertices_flag) { 1215 get_faces = PETSC_FALSE; 1216 get_edges = PETSC_FALSE; 1217 } 1218 if (pcbddc->constraints_flag) { 1219 get_vertices = PETSC_FALSE; 1220 } 1221 if (pcbddc->faces_flag) { 1222 get_edges = PETSC_FALSE; 1223 } 1224 if (pcbddc->edges_flag) { 1225 get_faces = PETSC_FALSE; 1226 } 1227 /* default */ 1228 if (!get_faces && !get_edges && !get_vertices) { 1229 get_vertices = PETSC_TRUE; 1230 } 1231 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,get_faces,get_edges,get_vertices,&n_ISForFaces,&ISForFaces,&n_ISForEdges,&ISForEdges,&ISForVertices); 1232 /* print some info */ 1233 if (pcbddc->dbg_flag) { 1234 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 1235 i = 0; 1236 if (ISForVertices) { 1237 ierr = ISGetSize(ISForVertices,&i);CHKERRQ(ierr); 1238 } 1239 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate vertices\n",PetscGlobalRank,i);CHKERRQ(ierr); 1240 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate edges\n",PetscGlobalRank,n_ISForEdges);CHKERRQ(ierr); 1241 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate faces\n",PetscGlobalRank,n_ISForFaces);CHKERRQ(ierr); 1242 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 1243 } 1244 /* check if near null space is attached to global mat */ 1245 ierr = MatGetNearNullSpace(pc->pmat,&nearnullsp);CHKERRQ(ierr); 1246 if (nearnullsp) { 1247 ierr = MatNullSpaceGetVecs(nearnullsp,&nnsp_has_cnst,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 1248 } else { /* if near null space is not provided BDDC uses constants by default */ 1249 nnsp_size = 0; 1250 nnsp_has_cnst = PETSC_TRUE; 1251 } 1252 /* get max number of constraints on a single cc */ 1253 max_constraints = nnsp_size; 1254 if (nnsp_has_cnst) max_constraints++; 1255 1256 /* 1257 Evaluate maximum storage size needed by the procedure 1258 - temp_indices will contain start index of each constraint stored as follows 1259 - temp_indices_to_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in local numbering) on which the constraint acts 1260 - temp_indices_to_constraint_B[temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in boundary numbering) on which the constraint acts 1261 - temp_quadrature_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the scalars representing the constraint itself 1262 */ 1263 total_counts = n_ISForFaces+n_ISForEdges; 1264 total_counts *= max_constraints; 1265 n_vertices = 0; 1266 if (ISForVertices) { 1267 ierr = ISGetSize(ISForVertices,&n_vertices);CHKERRQ(ierr); 1268 } 1269 total_counts += n_vertices; 1270 ierr = PetscMalloc((total_counts+1)*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 1271 ierr = PetscMalloc((total_counts+1)*sizeof(PetscBool),&change_basis);CHKERRQ(ierr); 1272 total_counts = 0; 1273 max_size_of_constraint = 0; 1274 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 1275 if (i<n_ISForEdges) { 1276 used_IS = &ISForEdges[i]; 1277 } else { 1278 used_IS = &ISForFaces[i-n_ISForEdges]; 1279 } 1280 ierr = ISGetSize(*used_IS,&j);CHKERRQ(ierr); 1281 total_counts += j; 1282 max_size_of_constraint = PetscMax(j,max_size_of_constraint); 1283 } 1284 total_counts *= max_constraints; 1285 total_counts += n_vertices; 1286 ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&temp_quadrature_constraint);CHKERRQ(ierr); 1287 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint);CHKERRQ(ierr); 1288 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint_B);CHKERRQ(ierr); 1289 /* local to boundary numbering */ 1290 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&local_to_B);CHKERRQ(ierr); 1291 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1292 for (i=0;i<pcis->n;i++) local_to_B[i]=-1; 1293 for (i=0;i<pcis->n_B;i++) local_to_B[is_indices[i]]=i; 1294 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1295 /* get local part of global near null space vectors */ 1296 ierr = PetscMalloc(nnsp_size*sizeof(Vec),&localnearnullsp);CHKERRQ(ierr); 1297 for (k=0;k<nnsp_size;k++) { 1298 ierr = VecDuplicate(pcis->vec1_N,&localnearnullsp[k]);CHKERRQ(ierr); 1299 ierr = VecScatterBegin(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1300 ierr = VecScatterEnd(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1301 } 1302 1303 /* whether or not to skip lapack calls */ 1304 skip_lapack = PETSC_TRUE; 1305 if (n_ISForFaces+n_ISForEdges) skip_lapack = PETSC_FALSE; 1306 1307 /* First we issue queries to allocate optimal workspace for LAPACKgesvd (or LAPACKsyev if SVD is missing) */ 1308 if (!pcbddc->use_nnsp_true && !skip_lapack) { 1309 PetscScalar temp_work; 1310 #if defined(PETSC_MISSING_LAPACK_GESVD) 1311 /* Proper Orthogonal Decomposition (POD) using the snapshot method */ 1312 ierr = PetscMalloc(max_constraints*max_constraints*sizeof(PetscScalar),&correlation_mat);CHKERRQ(ierr); 1313 ierr = PetscMalloc(max_constraints*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 1314 ierr = PetscMalloc(max_size_of_constraint*max_constraints*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr); 1315 #if defined(PETSC_USE_COMPLEX) 1316 ierr = PetscMalloc(3*max_constraints*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 1317 #endif 1318 /* now we evaluate the optimal workspace using query with lwork=-1 */ 1319 ierr = PetscBLASIntCast(max_constraints,&Blas_N);CHKERRQ(ierr); 1320 ierr = PetscBLASIntCast(max_constraints,&Blas_LDA);CHKERRQ(ierr); 1321 lwork = -1; 1322 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1323 #if !defined(PETSC_USE_COMPLEX) 1324 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,&temp_work,&lwork,&lierr)); 1325 #else 1326 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,&temp_work,&lwork,rwork,&lierr)); 1327 #endif 1328 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1329 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEV Lapack routine %d",(int)lierr); 1330 #else /* on missing GESVD */ 1331 /* SVD */ 1332 PetscInt max_n,min_n; 1333 max_n = max_size_of_constraint; 1334 min_n = max_constraints; 1335 if (max_size_of_constraint < max_constraints) { 1336 min_n = max_size_of_constraint; 1337 max_n = max_constraints; 1338 } 1339 ierr = PetscMalloc(min_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 1340 #if defined(PETSC_USE_COMPLEX) 1341 ierr = PetscMalloc(5*min_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 1342 #endif 1343 /* now we evaluate the optimal workspace using query with lwork=-1 */ 1344 lwork = -1; 1345 ierr = PetscBLASIntCast(max_n,&Blas_M);CHKERRQ(ierr); 1346 ierr = PetscBLASIntCast(min_n,&Blas_N);CHKERRQ(ierr); 1347 ierr = PetscBLASIntCast(max_n,&Blas_LDA);CHKERRQ(ierr); 1348 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1349 #if !defined(PETSC_USE_COMPLEX) 1350 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Blas_M,&Blas_N,&temp_quadrature_constraint[0],&Blas_LDA,singular_vals,&dummy_scalar_1,&dummy_int_1,&dummy_scalar_2,&dummy_int_2,&temp_work,&lwork,&lierr)); 1351 #else 1352 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Blas_M,&Blas_N,&temp_quadrature_constraint[0],&Blas_LDA,singular_vals,&dummy_scalar_1,&dummy_int_1,&dummy_scalar_2,&dummy_int_2,&temp_work,&lwork,rwork,&lierr)); 1353 #endif 1354 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1355 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GESVD Lapack routine %d",(int)lierr); 1356 #endif /* on missing GESVD */ 1357 /* Allocate optimal workspace */ 1358 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(temp_work),&lwork);CHKERRQ(ierr); 1359 ierr = PetscMalloc((PetscInt)lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr); 1360 } 1361 /* Now we can loop on constraining sets */ 1362 total_counts = 0; 1363 temp_indices[0] = 0; 1364 /* vertices */ 1365 if (ISForVertices) { 1366 ierr = ISGetIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1367 if (nnsp_has_cnst) { /* consider all vertices */ 1368 for (i=0;i<n_vertices;i++) { 1369 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 1370 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 1371 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 1372 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 1373 change_basis[total_counts]=PETSC_FALSE; 1374 total_counts++; 1375 } 1376 } else { /* consider vertices for which exist at least a localnearnullsp which is not null there */ 1377 PetscBool used_vertex; 1378 for (i=0;i<n_vertices;i++) { 1379 used_vertex = PETSC_FALSE; 1380 k = 0; 1381 while (!used_vertex && k<nnsp_size) { 1382 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 1383 if (PetscAbsScalar(array[is_indices[i]])>0.0) { 1384 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 1385 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 1386 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 1387 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 1388 change_basis[total_counts]=PETSC_FALSE; 1389 total_counts++; 1390 used_vertex = PETSC_TRUE; 1391 } 1392 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 1393 k++; 1394 } 1395 } 1396 } 1397 ierr = ISRestoreIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1398 n_vertices = total_counts; 1399 } 1400 1401 /* edges and faces */ 1402 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 1403 if (i<n_ISForEdges) { 1404 used_IS = &ISForEdges[i]; 1405 boolforchange = pcbddc->use_change_of_basis; /* change or not the basis on the edge */ 1406 } else { 1407 used_IS = &ISForFaces[i-n_ISForEdges]; 1408 boolforchange = (PetscBool)(pcbddc->use_change_of_basis && pcbddc->use_change_on_faces); /* change or not the basis on the face */ 1409 } 1410 temp_constraints = 0; /* zero the number of constraints I have on this conn comp */ 1411 temp_start_ptr = total_counts; /* need to know the starting index of constraints stored */ 1412 ierr = ISGetSize(*used_IS,&size_of_constraint);CHKERRQ(ierr); 1413 ierr = ISGetIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1414 /* change of basis should not be performed on local periodic nodes */ 1415 if (pcbddc->mat_graph->mirrors && pcbddc->mat_graph->mirrors[is_indices[0]]) boolforchange = PETSC_FALSE; 1416 if (nnsp_has_cnst) { 1417 PetscScalar quad_value; 1418 temp_constraints++; 1419 quad_value = (PetscScalar)(1.0/PetscSqrtReal((PetscReal)size_of_constraint)); 1420 for (j=0;j<size_of_constraint;j++) { 1421 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 1422 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 1423 temp_quadrature_constraint[temp_indices[total_counts]+j]=quad_value; 1424 } 1425 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 1426 change_basis[total_counts]=boolforchange; 1427 total_counts++; 1428 } 1429 for (k=0;k<nnsp_size;k++) { 1430 PetscReal real_value; 1431 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 1432 for (j=0;j<size_of_constraint;j++) { 1433 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 1434 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 1435 temp_quadrature_constraint[temp_indices[total_counts]+j]=array[is_indices[j]]; 1436 } 1437 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 1438 /* check if array is null on the connected component */ 1439 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1440 PetscStackCallBLAS("BLASasum",real_value = BLASasum_(&Blas_N,&temp_quadrature_constraint[temp_indices[total_counts]],&Blas_one)); 1441 if (real_value > 0.0) { /* keep indices and values */ 1442 temp_constraints++; 1443 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 1444 change_basis[total_counts]=boolforchange; 1445 total_counts++; 1446 } 1447 } 1448 ierr = ISRestoreIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 1449 /* perform SVD on the constraints if use_nnsp_true has not be requested by the user */ 1450 if (!pcbddc->use_nnsp_true) { 1451 PetscReal tol = 1.0e-8; /* tolerance for retaining eigenmodes */ 1452 1453 #if defined(PETSC_MISSING_LAPACK_GESVD) 1454 /* SVD: Y = U*S*V^H -> U (eigenvectors of Y*Y^H) = Y*V*(S)^\dag 1455 POD: Y^H*Y = V*D*V^H, D = S^H*S -> U = Y*V*D^(-1/2) 1456 -> When PETSC_USE_COMPLEX and PETSC_MISSING_LAPACK_GESVD are defined 1457 the constraints basis will differ (by a complex factor with absolute value equal to 1) 1458 from that computed using LAPACKgesvd 1459 -> This is due to a different computation of eigenvectors in LAPACKheev 1460 -> The quality of the POD-computed basis will be the same */ 1461 ierr = PetscMemzero(correlation_mat,temp_constraints*temp_constraints*sizeof(PetscScalar));CHKERRQ(ierr); 1462 /* Store upper triangular part of correlation matrix */ 1463 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1464 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1465 for (j=0;j<temp_constraints;j++) { 1466 for (k=0;k<j+1;k++) { 1467 PetscStackCallBLAS("BLASdot",correlation_mat[j*temp_constraints+k]=BLASdot_(&Blas_N,&temp_quadrature_constraint[temp_indices[temp_start_ptr+k]],&Blas_one,&temp_quadrature_constraint[temp_indices[temp_start_ptr+j]],&Blas_one_2)); 1468 } 1469 } 1470 /* compute eigenvalues and eigenvectors of correlation matrix */ 1471 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 1472 ierr = PetscBLASIntCast(temp_constraints,&Blas_LDA);CHKERRQ(ierr); 1473 #if !defined(PETSC_USE_COMPLEX) 1474 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,work,&lwork,&lierr)); 1475 #else 1476 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,work,&lwork,rwork,&lierr)); 1477 #endif 1478 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1479 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine %d",(int)lierr); 1480 /* retain eigenvalues greater than tol: note that LAPACKsyev gives eigs in ascending order */ 1481 j=0; 1482 while (j < temp_constraints && singular_vals[j] < tol) j++; 1483 total_counts=total_counts-j; 1484 /* scale and copy POD basis into used quadrature memory */ 1485 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1486 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 1487 ierr = PetscBLASIntCast(temp_constraints,&Blas_K);CHKERRQ(ierr); 1488 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1489 ierr = PetscBLASIntCast(temp_constraints,&Blas_LDB);CHKERRQ(ierr); 1490 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDC);CHKERRQ(ierr); 1491 if (j<temp_constraints) { 1492 PetscInt ii; 1493 for (k=j;k<temp_constraints;k++) singular_vals[k]=1.0/PetscSqrtReal(singular_vals[k]); 1494 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1495 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&Blas_M,&Blas_N,&Blas_K,&one,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Blas_LDA,correlation_mat,&Blas_LDB,&zero,temp_basis,&Blas_LDC)); 1496 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1497 for (k=0;k<temp_constraints-j;k++) { 1498 for (ii=0;ii<size_of_constraint;ii++) { 1499 temp_quadrature_constraint[temp_indices[temp_start_ptr+k]+ii]=singular_vals[temp_constraints-1-k]*temp_basis[(temp_constraints-1-k)*size_of_constraint+ii]; 1500 } 1501 } 1502 } 1503 #else /* on missing GESVD */ 1504 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1505 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 1506 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1507 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1508 #if !defined(PETSC_USE_COMPLEX) 1509 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Blas_M,&Blas_N,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Blas_LDA,singular_vals,&dummy_scalar_1,&dummy_int_1,&dummy_scalar_2,&dummy_int_2,work,&lwork,&lierr)); 1510 #else 1511 PetscStackCallBLAS("LAPACKgesvd",LAPACKgesvd_("O","N",&Blas_M,&Blas_N,&temp_quadrature_constraint[temp_indices[temp_start_ptr]],&Blas_LDA,singular_vals,&dummy_scalar_1,&dummy_int_1,&dummy_scalar_2,&dummy_int_2,work,&lwork,rwork,&lierr)); 1512 #endif 1513 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GESVD Lapack routine %d",(int)lierr); 1514 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1515 /* retain eigenvalues greater than tol: note that LAPACKgesvd gives eigs in descending order */ 1516 k = temp_constraints; 1517 if (k > size_of_constraint) k = size_of_constraint; 1518 j = 0; 1519 while (j < k && singular_vals[k-j-1] < tol) j++; 1520 total_counts = total_counts-temp_constraints+k-j; 1521 #endif /* on missing GESVD */ 1522 } 1523 } 1524 /* free index sets of faces, edges and vertices */ 1525 for (i=0;i<n_ISForFaces;i++) { 1526 ierr = ISDestroy(&ISForFaces[i]);CHKERRQ(ierr); 1527 } 1528 ierr = PetscFree(ISForFaces);CHKERRQ(ierr); 1529 for (i=0;i<n_ISForEdges;i++) { 1530 ierr = ISDestroy(&ISForEdges[i]);CHKERRQ(ierr); 1531 } 1532 ierr = PetscFree(ISForEdges);CHKERRQ(ierr); 1533 ierr = ISDestroy(&ISForVertices);CHKERRQ(ierr); 1534 1535 /* free workspace */ 1536 if (!pcbddc->use_nnsp_true && !skip_lapack) { 1537 ierr = PetscFree(work);CHKERRQ(ierr); 1538 #if defined(PETSC_USE_COMPLEX) 1539 ierr = PetscFree(rwork);CHKERRQ(ierr); 1540 #endif 1541 ierr = PetscFree(singular_vals);CHKERRQ(ierr); 1542 #if defined(PETSC_MISSING_LAPACK_GESVD) 1543 ierr = PetscFree(correlation_mat);CHKERRQ(ierr); 1544 ierr = PetscFree(temp_basis);CHKERRQ(ierr); 1545 #endif 1546 } 1547 for (k=0;k<nnsp_size;k++) { 1548 ierr = VecDestroy(&localnearnullsp[k]);CHKERRQ(ierr); 1549 } 1550 ierr = PetscFree(localnearnullsp);CHKERRQ(ierr); 1551 1552 /* set quantities in pcbddc data structure */ 1553 /* n_vertices defines the number of subdomain corners in the primal space */ 1554 /* n_constraints defines the number of averages (they can be point primal dofs if change of basis is requested) */ 1555 pcbddc->local_primal_size = total_counts; 1556 pcbddc->n_vertices = n_vertices; 1557 pcbddc->n_constraints = pcbddc->local_primal_size-pcbddc->n_vertices; 1558 1559 /* Create constraint matrix */ 1560 /* The constraint matrix is used to compute the l2g map of primal dofs */ 1561 /* so we need to set it up properly either with or without change of basis */ 1562 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ConstraintMatrix);CHKERRQ(ierr); 1563 ierr = MatSetType(pcbddc->ConstraintMatrix,impMatType);CHKERRQ(ierr); 1564 ierr = MatSetSizes(pcbddc->ConstraintMatrix,pcbddc->local_primal_size,pcis->n,pcbddc->local_primal_size,pcis->n);CHKERRQ(ierr); 1565 /* array to compute a local numbering of constraints : vertices first then constraints */ 1566 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&aux_primal_numbering);CHKERRQ(ierr); 1567 /* array to select the proper local node (of minimum index with respect to global ordering) when changing the basis */ 1568 /* note: it should not be needed since IS for faces and edges are already sorted by global ordering when analyzing the graph but... just in case */ 1569 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&aux_primal_minloc);CHKERRQ(ierr); 1570 /* auxiliary stuff for basis change */ 1571 ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&global_indices);CHKERRQ(ierr); 1572 ierr = PetscMalloc(pcis->n_B*sizeof(PetscBool),&touched);CHKERRQ(ierr); 1573 ierr = PetscMemzero(touched,pcis->n_B*sizeof(PetscBool));CHKERRQ(ierr); 1574 1575 /* find primal_dofs: subdomain corners plus dofs selected as primal after change of basis */ 1576 total_primal_vertices=0; 1577 for (i=0;i<pcbddc->local_primal_size;i++) { 1578 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 1579 if (size_of_constraint == 1) { 1580 touched[temp_indices_to_constraint_B[temp_indices[i]]]=PETSC_TRUE; 1581 aux_primal_numbering[total_primal_vertices]=temp_indices_to_constraint[temp_indices[i]]; 1582 aux_primal_minloc[total_primal_vertices]=0; 1583 total_primal_vertices++; 1584 } else if (change_basis[i]) { /* Same procedure used in PCBDDCGetPrimalConstraintsLocalIdx */ 1585 PetscInt min_loc,min_index; 1586 ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,&temp_indices_to_constraint[temp_indices[i]],global_indices);CHKERRQ(ierr); 1587 /* find first untouched local node */ 1588 k = 0; 1589 while (touched[temp_indices_to_constraint_B[temp_indices[i]+k]]) k++; 1590 min_index = global_indices[k]; 1591 min_loc = k; 1592 /* search the minimum among global nodes already untouched on the cc */ 1593 for (k=1;k<size_of_constraint;k++) { 1594 /* there can be more than one constraint on a single connected component */ 1595 if (min_index > global_indices[k] && !touched[temp_indices_to_constraint_B[temp_indices[i]+k]]) { 1596 min_index = global_indices[k]; 1597 min_loc = k; 1598 } 1599 } 1600 touched[temp_indices_to_constraint_B[temp_indices[i]+min_loc]] = PETSC_TRUE; 1601 aux_primal_numbering[total_primal_vertices]=temp_indices_to_constraint[temp_indices[i]+min_loc]; 1602 aux_primal_minloc[total_primal_vertices]=min_loc; 1603 total_primal_vertices++; 1604 } 1605 } 1606 /* free workspace */ 1607 ierr = PetscFree(global_indices);CHKERRQ(ierr); 1608 ierr = PetscFree(touched);CHKERRQ(ierr); 1609 /* permute indices in order to have a sorted set of vertices */ 1610 ierr = PetscSortInt(total_primal_vertices,aux_primal_numbering); 1611 1612 /* nonzero structure of constraint matrix */ 1613 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1614 for (i=0;i<total_primal_vertices;i++) nnz[i]=1; 1615 j=total_primal_vertices; 1616 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1617 if (!change_basis[i]) { 1618 nnz[j]=temp_indices[i+1]-temp_indices[i]; 1619 j++; 1620 } 1621 } 1622 ierr = MatSeqAIJSetPreallocation(pcbddc->ConstraintMatrix,0,nnz);CHKERRQ(ierr); 1623 ierr = PetscFree(nnz);CHKERRQ(ierr); 1624 /* set values in constraint matrix */ 1625 for (i=0;i<total_primal_vertices;i++) { 1626 ierr = MatSetValue(pcbddc->ConstraintMatrix,i,aux_primal_numbering[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 1627 } 1628 total_counts = total_primal_vertices; 1629 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1630 if (!change_basis[i]) { 1631 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 1632 ierr = MatSetValues(pcbddc->ConstraintMatrix,1,&total_counts,size_of_constraint,&temp_indices_to_constraint[temp_indices[i]],&temp_quadrature_constraint[temp_indices[i]],INSERT_VALUES);CHKERRQ(ierr); 1633 total_counts++; 1634 } 1635 } 1636 /* assembling */ 1637 ierr = MatAssemblyBegin(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1638 ierr = MatAssemblyEnd(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1639 /* 1640 ierr = MatView(pcbddc->ConstraintMatrix,(PetscViewer)0);CHKERRQ(ierr); 1641 */ 1642 /* Create matrix for change of basis. We don't need it in case pcbddc->use_change_of_basis is FALSE */ 1643 if (pcbddc->use_change_of_basis) { 1644 PetscBool qr_needed = PETSC_FALSE; 1645 /* change of basis acts on local interfaces -> dimension is n_B x n_B */ 1646 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1647 ierr = MatSetType(pcbddc->ChangeOfBasisMatrix,impMatType);CHKERRQ(ierr); 1648 ierr = MatSetSizes(pcbddc->ChangeOfBasisMatrix,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B);CHKERRQ(ierr); 1649 /* work arrays */ 1650 ierr = PetscMalloc(pcis->n_B*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1651 for (i=0;i<pcis->n_B;i++) nnz[i]=1; 1652 /* nonzeros per row */ 1653 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1654 if (change_basis[i]) { 1655 qr_needed = PETSC_TRUE; 1656 size_of_constraint = temp_indices[i+1]-temp_indices[i]; 1657 for (j=0;j<size_of_constraint;j++) nnz[temp_indices_to_constraint_B[temp_indices[i]+j]] = size_of_constraint; 1658 } 1659 } 1660 ierr = MatSeqAIJSetPreallocation(pcbddc->ChangeOfBasisMatrix,0,nnz);CHKERRQ(ierr); 1661 ierr = PetscFree(nnz);CHKERRQ(ierr); 1662 /* Set initial identity in the matrix */ 1663 for (i=0;i<pcis->n_B;i++) { 1664 ierr = MatSetValue(pcbddc->ChangeOfBasisMatrix,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr); 1665 } 1666 1667 /* Now we loop on the constraints which need a change of basis */ 1668 /* Change of basis matrix is evaluated as the FIRST APPROACH in */ 1669 /* Klawonn and Widlund, Dual-primal FETI-DP methods for linear elasticity, (see Sect 6.2.1) */ 1670 /* Change of basis matrix T computed via QR decomposition of constraints */ 1671 if (qr_needed) { 1672 /* dual and primal dofs on a single cc */ 1673 PetscInt dual_dofs,primal_dofs; 1674 /* iterator on aux_primal_minloc (ordered as read from nearnullspace: vertices, edges and then constraints) */ 1675 PetscInt primal_counter; 1676 /* working stuff for GEQRF */ 1677 PetscScalar *qr_basis,*qr_tau,*qr_work,lqr_work_t; 1678 PetscBLASInt lqr_work; 1679 /* working stuff for UNGQR */ 1680 PetscScalar *gqr_work,lgqr_work_t; 1681 PetscBLASInt lgqr_work; 1682 /* working stuff for TRTRS */ 1683 PetscScalar *trs_rhs; 1684 PetscBLASInt Blas_NRHS; 1685 /* pointers for values insertion into change of basis matrix */ 1686 PetscInt *start_rows,*start_cols; 1687 PetscScalar *start_vals; 1688 /* working stuff for values insertion */ 1689 PetscBool *is_primal; 1690 1691 /* space to store Q */ 1692 ierr = PetscMalloc((max_size_of_constraint)*(max_size_of_constraint)*sizeof(PetscScalar),&qr_basis);CHKERRQ(ierr); 1693 /* first we issue queries for optimal work */ 1694 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_M);CHKERRQ(ierr); 1695 ierr = PetscBLASIntCast(max_constraints,&Blas_N);CHKERRQ(ierr); 1696 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1697 lqr_work = -1; 1698 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,&lqr_work_t,&lqr_work,&lierr)); 1699 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GEQRF Lapack routine %d",(int)lierr); 1700 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lqr_work_t),&lqr_work);CHKERRQ(ierr); 1701 ierr = PetscMalloc((PetscInt)PetscRealPart(lqr_work_t)*sizeof(*qr_work),&qr_work);CHKERRQ(ierr); 1702 lgqr_work = -1; 1703 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_M);CHKERRQ(ierr); 1704 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_N);CHKERRQ(ierr); 1705 ierr = PetscBLASIntCast(max_constraints,&Blas_K);CHKERRQ(ierr); 1706 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1707 if (Blas_K>Blas_M) Blas_K=Blas_M; /* adjust just for computing optimal work */ 1708 PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Blas_M,&Blas_N,&Blas_K,qr_basis,&Blas_LDA,qr_tau,&lgqr_work_t,&lgqr_work,&lierr)); 1709 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to UNGQR Lapack routine %d",(int)lierr); 1710 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lgqr_work_t),&lgqr_work);CHKERRQ(ierr); 1711 ierr = PetscMalloc((PetscInt)PetscRealPart(lgqr_work_t)*sizeof(*gqr_work),&gqr_work);CHKERRQ(ierr); 1712 /* array to store scaling factors for reflectors */ 1713 ierr = PetscMalloc(max_constraints*sizeof(*qr_tau),&qr_tau);CHKERRQ(ierr); 1714 /* array to store rhs and solution of triangular solver */ 1715 ierr = PetscMalloc(max_constraints*max_constraints*sizeof(*trs_rhs),&trs_rhs);CHKERRQ(ierr); 1716 /* array to store whether a node is primal or not */ 1717 ierr = PetscMalloc(pcis->n_B*sizeof(*is_primal),&is_primal);CHKERRQ(ierr); 1718 ierr = PetscMemzero(is_primal,pcis->n_B*sizeof(*is_primal));CHKERRQ(ierr); 1719 for (i=0;i<total_primal_vertices;i++) is_primal[local_to_B[aux_primal_numbering[i]]] = PETSC_TRUE; 1720 1721 /* allocating workspace for check */ 1722 if (pcbddc->dbg_flag) { 1723 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 1724 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Checking change of basis computation for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 1725 ierr = PetscMalloc(max_size_of_constraint*(max_constraints+max_size_of_constraint)*sizeof(*work),&work);CHKERRQ(ierr); 1726 } 1727 1728 /* loop on constraints and see whether or not they need a change of basis */ 1729 /* -> using implicit ordering contained in temp_indices data */ 1730 total_counts = pcbddc->n_vertices; 1731 primal_counter = total_counts; 1732 while (total_counts<pcbddc->local_primal_size) { 1733 primal_dofs = 1; 1734 if (change_basis[total_counts]) { 1735 /* get all constraints with same support: if more then one constraint is present on the cc then surely indices are stored contiguosly */ 1736 while (total_counts+primal_dofs < pcbddc->local_primal_size && temp_indices_to_constraint_B[temp_indices[total_counts]] == temp_indices_to_constraint_B[temp_indices[total_counts+primal_dofs]]) { 1737 primal_dofs++; 1738 } 1739 /* get constraint info */ 1740 size_of_constraint = temp_indices[total_counts+1]-temp_indices[total_counts]; 1741 dual_dofs = size_of_constraint-primal_dofs; 1742 1743 /* copy quadrature constraints for change of basis check */ 1744 if (pcbddc->dbg_flag) { 1745 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Constraint %d to %d need a change of basis (size %d)\n",total_counts,total_counts+primal_dofs,size_of_constraint);CHKERRQ(ierr); 1746 ierr = PetscMemcpy(work,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1747 } 1748 1749 /* copy temporary constraints into larger work vector (in order to store all columns of Q) */ 1750 ierr = PetscMemcpy(qr_basis,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1751 1752 /* compute QR decomposition of constraints */ 1753 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1754 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1755 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1756 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1757 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,qr_work,&lqr_work,&lierr)); 1758 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GEQRF Lapack routine %d",(int)lierr); 1759 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1760 1761 /* explictly compute R^-T */ 1762 ierr = PetscMemzero(trs_rhs,primal_dofs*primal_dofs*sizeof(*trs_rhs));CHKERRQ(ierr); 1763 for (j=0;j<primal_dofs;j++) trs_rhs[j*(primal_dofs+1)] = 1.0; 1764 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1765 ierr = PetscBLASIntCast(primal_dofs,&Blas_NRHS);CHKERRQ(ierr); 1766 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1767 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDB);CHKERRQ(ierr); 1768 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1769 PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&Blas_N,&Blas_NRHS,qr_basis,&Blas_LDA,trs_rhs,&Blas_LDB,&lierr)); 1770 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TRTRS Lapack routine %d",(int)lierr); 1771 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1772 1773 /* explcitly compute all columns of Q (Q = [Q1 | Q2] ) overwriting QR factorization in qr_basis */ 1774 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1775 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1776 ierr = PetscBLASIntCast(primal_dofs,&Blas_K);CHKERRQ(ierr); 1777 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1778 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1779 PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Blas_M,&Blas_N,&Blas_K,qr_basis,&Blas_LDA,qr_tau,gqr_work,&lgqr_work,&lierr)); 1780 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in UNGQR Lapack routine %d",(int)lierr); 1781 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1782 1783 /* first primal_dofs columns of Q need to be re-scaled in order to be unitary w.r.t constraints 1784 i.e. C_{pxn}*Q_{nxn} should be equal to [I_pxp | 0_pxd] (see check below) 1785 where n=size_of_constraint, p=primal_dofs, d=dual_dofs (n=p+d), I and 0 identity and null matrix resp. */ 1786 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1787 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1788 ierr = PetscBLASIntCast(primal_dofs,&Blas_K);CHKERRQ(ierr); 1789 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1790 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDB);CHKERRQ(ierr); 1791 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDC);CHKERRQ(ierr); 1792 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1793 PetscStackCallBLAS("BLASgemm",BLASgemm_("N","N",&Blas_M,&Blas_N,&Blas_K,&one,qr_basis,&Blas_LDA,trs_rhs,&Blas_LDB,&zero,&temp_quadrature_constraint[temp_indices[total_counts]],&Blas_LDC)); 1794 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1795 ierr = PetscMemcpy(qr_basis,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1796 1797 /* insert values in change of basis matrix respecting global ordering of new primal dofs */ 1798 start_rows = &temp_indices_to_constraint_B[temp_indices[total_counts]]; 1799 /* insert cols for primal dofs */ 1800 for (j=0;j<primal_dofs;j++) { 1801 start_vals = &qr_basis[j*size_of_constraint]; 1802 start_cols = &temp_indices_to_constraint_B[temp_indices[total_counts]+aux_primal_minloc[primal_counter+j]]; 1803 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,start_rows,1,start_cols,start_vals,INSERT_VALUES);CHKERRQ(ierr); 1804 } 1805 /* insert cols for dual dofs */ 1806 for (j=0,k=0;j<dual_dofs;k++) { 1807 if (!is_primal[temp_indices_to_constraint_B[temp_indices[total_counts]+k]]) { 1808 start_vals = &qr_basis[(primal_dofs+j)*size_of_constraint]; 1809 start_cols = &temp_indices_to_constraint_B[temp_indices[total_counts]+k]; 1810 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,start_rows,1,start_cols,start_vals,INSERT_VALUES);CHKERRQ(ierr); 1811 j++; 1812 } 1813 } 1814 1815 /* check change of basis */ 1816 if (pcbddc->dbg_flag) { 1817 PetscInt ii,jj; 1818 PetscBool valid_qr=PETSC_TRUE; 1819 ierr = PetscBLASIntCast(primal_dofs,&Blas_M);CHKERRQ(ierr); 1820 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1821 ierr = PetscBLASIntCast(size_of_constraint,&Blas_K);CHKERRQ(ierr); 1822 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1823 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDB);CHKERRQ(ierr); 1824 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDC);CHKERRQ(ierr); 1825 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1826 PetscStackCallBLAS("BLASgemm",BLASgemm_("T","N",&Blas_M,&Blas_N,&Blas_K,&one,work,&Blas_LDA,qr_basis,&Blas_LDB,&zero,&work[size_of_constraint*primal_dofs],&Blas_LDC)); 1827 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1828 for (jj=0;jj<size_of_constraint;jj++) { 1829 for (ii=0;ii<primal_dofs;ii++) { 1830 if (ii != jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]) > 1.e-12) valid_qr = PETSC_FALSE; 1831 if (ii == jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]-1.0) > 1.e-12) valid_qr = PETSC_FALSE; 1832 } 1833 } 1834 if (!valid_qr) { 1835 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\t-> wrong change of basis!\n",PetscGlobalRank);CHKERRQ(ierr); 1836 for (jj=0;jj<size_of_constraint;jj++) { 1837 for (ii=0;ii<primal_dofs;ii++) { 1838 if (ii != jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]) > 1.e-12) { 1839 PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\tQr basis function %d is not orthogonal to constraint %d (%1.14e)!\n",jj,ii,PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii])); 1840 } 1841 if (ii == jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]-1.0) > 1.e-12) { 1842 PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\tQr basis function %d is not unitary w.r.t constraint %d (%1.14e)!\n",jj,ii,PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii])); 1843 } 1844 } 1845 } 1846 } else { 1847 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\t-> right change of basis!\n",PetscGlobalRank);CHKERRQ(ierr); 1848 } 1849 } 1850 /* increment primal counter */ 1851 primal_counter += primal_dofs; 1852 } else { 1853 if (pcbddc->dbg_flag) { 1854 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Constraint %d does not need a change of basis (size %d)\n",total_counts,temp_indices[total_counts+1]-temp_indices[total_counts]);CHKERRQ(ierr); 1855 } 1856 } 1857 /* increment constraint counter total_counts */ 1858 total_counts += primal_dofs; 1859 } 1860 if (pcbddc->dbg_flag) { 1861 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 1862 ierr = PetscFree(work);CHKERRQ(ierr); 1863 } 1864 /* free workspace */ 1865 ierr = PetscFree(trs_rhs);CHKERRQ(ierr); 1866 ierr = PetscFree(qr_tau);CHKERRQ(ierr); 1867 ierr = PetscFree(qr_work);CHKERRQ(ierr); 1868 ierr = PetscFree(gqr_work);CHKERRQ(ierr); 1869 ierr = PetscFree(is_primal);CHKERRQ(ierr); 1870 ierr = PetscFree(qr_basis);CHKERRQ(ierr); 1871 } 1872 /* assembling */ 1873 ierr = MatAssemblyBegin(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1874 ierr = MatAssemblyEnd(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1875 /* 1876 ierr = MatView(pcbddc->ChangeOfBasisMatrix,(PetscViewer)0);CHKERRQ(ierr); 1877 */ 1878 } 1879 /* free workspace */ 1880 ierr = PetscFree(aux_primal_numbering);CHKERRQ(ierr); 1881 ierr = PetscFree(aux_primal_minloc);CHKERRQ(ierr); 1882 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 1883 ierr = PetscFree(change_basis);CHKERRQ(ierr); 1884 ierr = PetscFree(temp_indices_to_constraint);CHKERRQ(ierr); 1885 ierr = PetscFree(temp_indices_to_constraint_B);CHKERRQ(ierr); 1886 ierr = PetscFree(local_to_B);CHKERRQ(ierr); 1887 ierr = PetscFree(temp_quadrature_constraint);CHKERRQ(ierr); 1888 PetscFunctionReturn(0); 1889 } 1890 1891 #undef __FUNCT__ 1892 #define __FUNCT__ "PCBDDCAnalyzeInterface" 1893 PetscErrorCode PCBDDCAnalyzeInterface(PC pc) 1894 { 1895 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1896 PC_IS *pcis = (PC_IS*)pc->data; 1897 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1898 PetscInt bs,ierr,i,vertex_size; 1899 PetscViewer viewer=pcbddc->dbg_viewer; 1900 1901 PetscFunctionBegin; 1902 /* Init local Graph struct */ 1903 ierr = PCBDDCGraphInit(pcbddc->mat_graph,matis->mapping);CHKERRQ(ierr); 1904 1905 /* Check validity of the csr graph passed in by the user */ 1906 if (pcbddc->mat_graph->nvtxs_csr != pcbddc->mat_graph->nvtxs) { 1907 ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr); 1908 } 1909 /* Set default CSR adjacency of local dofs if not provided by the user with PCBDDCSetLocalAdjacencyGraph */ 1910 if (!pcbddc->mat_graph->xadj || !pcbddc->mat_graph->adjncy) { 1911 Mat mat_adj; 1912 const PetscInt *xadj,*adjncy; 1913 PetscBool flg_row=PETSC_TRUE; 1914 1915 ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr); 1916 ierr = MatGetRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1917 if (!flg_row) { 1918 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called in %s\n",__FUNCT__); 1919 } 1920 ierr = PCBDDCSetLocalAdjacencyGraph(pc,i,xadj,adjncy,PETSC_COPY_VALUES);CHKERRQ(ierr); 1921 ierr = MatRestoreRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1922 if (!flg_row) { 1923 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__); 1924 } 1925 ierr = MatDestroy(&mat_adj);CHKERRQ(ierr); 1926 } 1927 1928 /* Set default dofs' splitting if no information has been provided by the user with PCBDDCSetDofsSplitting */ 1929 vertex_size = 1; 1930 if (!pcbddc->n_ISForDofs) { 1931 IS *custom_ISForDofs; 1932 1933 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1934 ierr = PetscMalloc(bs*sizeof(IS),&custom_ISForDofs);CHKERRQ(ierr); 1935 for (i=0;i<bs;i++) { 1936 ierr = ISCreateStride(PETSC_COMM_SELF,pcis->n/bs,i,bs,&custom_ISForDofs[i]);CHKERRQ(ierr); 1937 } 1938 ierr = PCBDDCSetDofsSplitting(pc,bs,custom_ISForDofs);CHKERRQ(ierr); 1939 /* remove my references to IS objects */ 1940 for (i=0;i<bs;i++) { 1941 ierr = ISDestroy(&custom_ISForDofs[i]);CHKERRQ(ierr); 1942 } 1943 ierr = PetscFree(custom_ISForDofs);CHKERRQ(ierr); 1944 } else { /* mat block size as vertex size (used for elasticity) */ 1945 ierr = MatGetBlockSize(matis->A,&vertex_size);CHKERRQ(ierr); 1946 } 1947 1948 /* Setup of Graph */ 1949 ierr = PCBDDCGraphSetUp(pcbddc->mat_graph,vertex_size,pcbddc->NeumannBoundaries,pcbddc->DirichletBoundaries,pcbddc->n_ISForDofs,pcbddc->ISForDofs,pcbddc->user_primal_vertices); 1950 1951 /* Graph's connected components analysis */ 1952 ierr = PCBDDCGraphComputeConnectedComponents(pcbddc->mat_graph);CHKERRQ(ierr); 1953 1954 /* print some info to stdout */ 1955 if (pcbddc->dbg_flag) { 1956 ierr = PCBDDCGraphASCIIView(pcbddc->mat_graph,pcbddc->dbg_flag,viewer); 1957 } 1958 PetscFunctionReturn(0); 1959 } 1960 1961 #undef __FUNCT__ 1962 #define __FUNCT__ "PCBDDCGetPrimalVerticesLocalIdx" 1963 PetscErrorCode PCBDDCGetPrimalVerticesLocalIdx(PC pc, PetscInt *n_vertices, PetscInt *vertices_idx[]) 1964 { 1965 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1966 PetscInt *vertices,*row_cmat_indices,n,i,size_of_constraint,local_primal_size; 1967 PetscErrorCode ierr; 1968 1969 PetscFunctionBegin; 1970 n = 0; 1971 vertices = 0; 1972 if (pcbddc->ConstraintMatrix) { 1973 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&i);CHKERRQ(ierr); 1974 for (i=0;i<local_primal_size;i++) { 1975 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1976 if (size_of_constraint == 1) n++; 1977 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1978 } 1979 if (vertices_idx) { 1980 ierr = PetscMalloc(n*sizeof(PetscInt),&vertices);CHKERRQ(ierr); 1981 n = 0; 1982 for (i=0;i<local_primal_size;i++) { 1983 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1984 if (size_of_constraint == 1) { 1985 vertices[n++]=row_cmat_indices[0]; 1986 } 1987 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1988 } 1989 } 1990 } 1991 *n_vertices = n; 1992 if (vertices_idx) *vertices_idx = vertices; 1993 PetscFunctionReturn(0); 1994 } 1995 1996 #undef __FUNCT__ 1997 #define __FUNCT__ "PCBDDCGetPrimalConstraintsLocalIdx" 1998 PetscErrorCode PCBDDCGetPrimalConstraintsLocalIdx(PC pc, PetscInt *n_constraints, PetscInt *constraints_idx[]) 1999 { 2000 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 2001 PetscInt *constraints_index,*row_cmat_indices,*row_cmat_global_indices; 2002 PetscInt n,i,j,size_of_constraint,local_primal_size,local_size,max_size_of_constraint,min_index,min_loc; 2003 PetscBool *touched; 2004 PetscErrorCode ierr; 2005 2006 PetscFunctionBegin; 2007 n = 0; 2008 constraints_index = 0; 2009 if (pcbddc->ConstraintMatrix) { 2010 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&local_size);CHKERRQ(ierr); 2011 max_size_of_constraint = 0; 2012 for (i=0;i<local_primal_size;i++) { 2013 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 2014 if (size_of_constraint > 1) { 2015 n++; 2016 } 2017 max_size_of_constraint = PetscMax(size_of_constraint,max_size_of_constraint); 2018 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 2019 } 2020 if (constraints_idx) { 2021 ierr = PetscMalloc(n*sizeof(PetscInt),&constraints_index);CHKERRQ(ierr); 2022 ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&row_cmat_global_indices);CHKERRQ(ierr); 2023 ierr = PetscMalloc(local_size*sizeof(PetscBool),&touched);CHKERRQ(ierr); 2024 ierr = PetscMemzero(touched,local_size*sizeof(PetscBool));CHKERRQ(ierr); 2025 n = 0; 2026 for (i=0;i<local_primal_size;i++) { 2027 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 2028 if (size_of_constraint > 1) { 2029 ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,row_cmat_indices,row_cmat_global_indices);CHKERRQ(ierr); 2030 /* find first untouched local node */ 2031 j = 0; 2032 while(touched[row_cmat_indices[j]]) j++; 2033 min_index = row_cmat_global_indices[j]; 2034 min_loc = j; 2035 /* search the minimum among nodes not yet touched on the connected component 2036 since there can be more than one constraint on a single cc */ 2037 for (j=1;j<size_of_constraint;j++) { 2038 if (min_index > row_cmat_global_indices[j] && !touched[row_cmat_indices[j]]) { 2039 min_index = row_cmat_global_indices[j]; 2040 min_loc = j; 2041 } 2042 } 2043 touched[row_cmat_indices[min_loc]] = PETSC_TRUE; 2044 constraints_index[n++] = row_cmat_indices[min_loc]; 2045 } 2046 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 2047 } 2048 ierr = PetscFree(touched);CHKERRQ(ierr); 2049 ierr = PetscFree(row_cmat_global_indices);CHKERRQ(ierr); 2050 } 2051 } 2052 *n_constraints = n; 2053 if (constraints_idx) *constraints_idx = constraints_index; 2054 PetscFunctionReturn(0); 2055 } 2056 2057 /* the next two functions has been adapted from pcis.c */ 2058 #undef __FUNCT__ 2059 #define __FUNCT__ "PCBDDCApplySchur" 2060 PetscErrorCode PCBDDCApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 2061 { 2062 PetscErrorCode ierr; 2063 PC_IS *pcis = (PC_IS*)(pc->data); 2064 2065 PetscFunctionBegin; 2066 if (!vec2_B) { vec2_B = v; } 2067 ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 2068 ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr); 2069 ierr = KSPSolve(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 2070 ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr); 2071 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 2072 PetscFunctionReturn(0); 2073 } 2074 2075 #undef __FUNCT__ 2076 #define __FUNCT__ "PCBDDCApplySchurTranspose" 2077 PetscErrorCode PCBDDCApplySchurTranspose(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 2078 { 2079 PetscErrorCode ierr; 2080 PC_IS *pcis = (PC_IS*)(pc->data); 2081 2082 PetscFunctionBegin; 2083 if (!vec2_B) { vec2_B = v; } 2084 ierr = MatMultTranspose(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 2085 ierr = MatMultTranspose(pcis->A_BI,v,vec1_D);CHKERRQ(ierr); 2086 ierr = KSPSolveTranspose(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 2087 ierr = MatMultTranspose(pcis->A_IB,vec2_D,vec2_B);CHKERRQ(ierr); 2088 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 2089 PetscFunctionReturn(0); 2090 } 2091 2092 #undef __FUNCT__ 2093 #define __FUNCT__ "PCBDDCSubsetNumbering" 2094 PetscErrorCode PCBDDCSubsetNumbering(MPI_Comm comm,ISLocalToGlobalMapping l2gmap, PetscInt n_local_dofs, PetscInt local_dofs[], PetscInt local_dofs_mult[], PetscInt* n_global_subset, PetscInt* global_numbering_subset[]) 2095 { 2096 Vec local_vec,global_vec; 2097 IS seqis,paris; 2098 VecScatter scatter_ctx; 2099 PetscScalar *array; 2100 PetscInt *temp_global_dofs; 2101 PetscScalar globalsum; 2102 PetscInt i,j,s; 2103 PetscInt nlocals,first_index,old_index,max_local; 2104 PetscMPIInt rank_prec_comm,size_prec_comm,max_global; 2105 PetscMPIInt *dof_sizes,*dof_displs; 2106 PetscBool first_found; 2107 PetscErrorCode ierr; 2108 2109 PetscFunctionBegin; 2110 /* mpi buffers */ 2111 MPI_Comm_size(comm,&size_prec_comm); 2112 MPI_Comm_rank(comm,&rank_prec_comm); 2113 j = ( !rank_prec_comm ? size_prec_comm : 0); 2114 ierr = PetscMalloc(j*sizeof(*dof_sizes),&dof_sizes);CHKERRQ(ierr); 2115 ierr = PetscMalloc(j*sizeof(*dof_displs),&dof_displs);CHKERRQ(ierr); 2116 /* get maximum size of subset */ 2117 ierr = PetscMalloc(n_local_dofs*sizeof(PetscInt),&temp_global_dofs);CHKERRQ(ierr); 2118 ierr = ISLocalToGlobalMappingApply(l2gmap,n_local_dofs,local_dofs,temp_global_dofs);CHKERRQ(ierr); 2119 max_local = 0; 2120 if (n_local_dofs) { 2121 max_local = temp_global_dofs[0]; 2122 for (i=1;i<n_local_dofs;i++) { 2123 if (max_local < temp_global_dofs[i] ) { 2124 max_local = temp_global_dofs[i]; 2125 } 2126 } 2127 } 2128 ierr = MPI_Allreduce(&max_local,&max_global,1,MPIU_INT,MPI_MAX,comm); 2129 max_global++; 2130 max_local = 0; 2131 if (n_local_dofs) { 2132 max_local = local_dofs[0]; 2133 for (i=1;i<n_local_dofs;i++) { 2134 if (max_local < local_dofs[i] ) { 2135 max_local = local_dofs[i]; 2136 } 2137 } 2138 } 2139 max_local++; 2140 /* allocate workspace */ 2141 ierr = VecCreate(PETSC_COMM_SELF,&local_vec);CHKERRQ(ierr); 2142 ierr = VecSetSizes(local_vec,PETSC_DECIDE,max_local);CHKERRQ(ierr); 2143 ierr = VecSetType(local_vec,VECSEQ);CHKERRQ(ierr); 2144 ierr = VecCreate(comm,&global_vec);CHKERRQ(ierr); 2145 ierr = VecSetSizes(global_vec,PETSC_DECIDE,max_global);CHKERRQ(ierr); 2146 ierr = VecSetType(global_vec,VECMPI);CHKERRQ(ierr); 2147 /* create scatter */ 2148 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_local_dofs,local_dofs,PETSC_COPY_VALUES,&seqis);CHKERRQ(ierr); 2149 ierr = ISCreateGeneral(comm,n_local_dofs,temp_global_dofs,PETSC_COPY_VALUES,&paris);CHKERRQ(ierr); 2150 ierr = VecScatterCreate(local_vec,seqis,global_vec,paris,&scatter_ctx);CHKERRQ(ierr); 2151 ierr = ISDestroy(&seqis);CHKERRQ(ierr); 2152 ierr = ISDestroy(&paris);CHKERRQ(ierr); 2153 /* init array */ 2154 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 2155 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 2156 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 2157 if (local_dofs_mult) { 2158 for (i=0;i<n_local_dofs;i++) { 2159 array[local_dofs[i]]=(PetscScalar)local_dofs_mult[i]; 2160 } 2161 } else { 2162 for (i=0;i<n_local_dofs;i++) { 2163 array[local_dofs[i]]=1.0; 2164 } 2165 } 2166 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 2167 /* scatter into global vec and get total number of global dofs */ 2168 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2169 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2170 ierr = VecSum(global_vec,&globalsum);CHKERRQ(ierr); 2171 *n_global_subset = (PetscInt)PetscRealPart(globalsum); 2172 /* Fill global_vec with cumulative function for global numbering */ 2173 ierr = VecGetArray(global_vec,&array);CHKERRQ(ierr); 2174 ierr = VecGetLocalSize(global_vec,&s);CHKERRQ(ierr); 2175 nlocals = 0; 2176 first_index = -1; 2177 first_found = PETSC_FALSE; 2178 for (i=0;i<s;i++) { 2179 if (!first_found && PetscRealPart(array[i]) > 0.0) { 2180 first_found = PETSC_TRUE; 2181 first_index = i; 2182 } 2183 nlocals += (PetscInt)PetscRealPart(array[i]); 2184 } 2185 ierr = MPI_Gather(&nlocals,1,MPIU_INT,dof_sizes,1,MPIU_INT,0,comm);CHKERRQ(ierr); 2186 if (!rank_prec_comm) { 2187 dof_displs[0]=0; 2188 for (i=1;i<size_prec_comm;i++) { 2189 dof_displs[i] = dof_displs[i-1]+dof_sizes[i-1]; 2190 } 2191 } 2192 ierr = MPI_Scatter(dof_displs,1,MPIU_INT,&nlocals,1,MPIU_INT,0,comm);CHKERRQ(ierr); 2193 if (first_found) { 2194 array[first_index] += (PetscScalar)nlocals; 2195 old_index = first_index; 2196 for (i=first_index+1;i<s;i++) { 2197 if (PetscRealPart(array[i]) > 0.0) { 2198 array[i] += array[old_index]; 2199 old_index = i; 2200 } 2201 } 2202 } 2203 ierr = VecRestoreArray(global_vec,&array);CHKERRQ(ierr); 2204 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 2205 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2206 ierr = VecScatterEnd (scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 2207 /* get global ordering of local dofs */ 2208 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 2209 if (local_dofs_mult) { 2210 for (i=0;i<n_local_dofs;i++) { 2211 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-local_dofs_mult[i]; 2212 } 2213 } else { 2214 for (i=0;i<n_local_dofs;i++) { 2215 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-1; 2216 } 2217 } 2218 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 2219 /* free workspace */ 2220 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 2221 ierr = VecDestroy(&local_vec);CHKERRQ(ierr); 2222 ierr = VecDestroy(&global_vec);CHKERRQ(ierr); 2223 ierr = PetscFree(dof_sizes);CHKERRQ(ierr); 2224 ierr = PetscFree(dof_displs);CHKERRQ(ierr); 2225 /* return pointer to global ordering of local dofs */ 2226 *global_numbering_subset = temp_global_dofs; 2227 PetscFunctionReturn(0); 2228 } 2229 2230 #undef __FUNCT__ 2231 #define __FUNCT__ "PCBDDCOrthonormalizeVecs" 2232 PetscErrorCode PCBDDCOrthonormalizeVecs(PetscInt n, Vec vecs[]) 2233 { 2234 PetscInt i,j; 2235 PetscScalar *alphas; 2236 PetscErrorCode ierr; 2237 2238 PetscFunctionBegin; 2239 /* this implements stabilized Gram-Schmidt */ 2240 ierr = PetscMalloc(n*sizeof(PetscScalar),&alphas);CHKERRQ(ierr); 2241 for (i=0;i<n;i++) { 2242 ierr = VecNormalize(vecs[i],NULL);CHKERRQ(ierr); 2243 if (i<n) { ierr = VecMDot(vecs[i],n-i-1,&vecs[i+1],&alphas[i+1]);CHKERRQ(ierr); } 2244 for (j=i+1;j<n;j++) { ierr = VecAXPY(vecs[j],PetscConj(-alphas[j]),vecs[i]);CHKERRQ(ierr); } 2245 } 2246 ierr = PetscFree(alphas);CHKERRQ(ierr); 2247 PetscFunctionReturn(0); 2248 } 2249 2250 2251