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