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->temp_solution);CHKERRQ(ierr); 49 ierr = VecDestroy(&pcbddc->original_rhs);CHKERRQ(ierr); 50 ierr = VecDestroy(&pcbddc->coarse_vec);CHKERRQ(ierr); 51 ierr = VecDestroy(&pcbddc->coarse_rhs);CHKERRQ(ierr); 52 ierr = KSPDestroy(&pcbddc->coarse_ksp);CHKERRQ(ierr); 53 ierr = MatDestroy(&pcbddc->coarse_mat);CHKERRQ(ierr); 54 ierr = MatDestroy(&pcbddc->coarse_phi_B);CHKERRQ(ierr); 55 ierr = MatDestroy(&pcbddc->coarse_phi_D);CHKERRQ(ierr); 56 ierr = MatDestroy(&pcbddc->coarse_psi_B);CHKERRQ(ierr); 57 ierr = MatDestroy(&pcbddc->coarse_psi_D);CHKERRQ(ierr); 58 ierr = VecDestroy(&pcbddc->vec1_P);CHKERRQ(ierr); 59 ierr = VecDestroy(&pcbddc->vec1_C);CHKERRQ(ierr); 60 ierr = MatDestroy(&pcbddc->local_auxmat1);CHKERRQ(ierr); 61 ierr = MatDestroy(&pcbddc->local_auxmat2);CHKERRQ(ierr); 62 ierr = VecDestroy(&pcbddc->vec1_R);CHKERRQ(ierr); 63 ierr = VecDestroy(&pcbddc->vec2_R);CHKERRQ(ierr); 64 ierr = VecDestroy(&pcbddc->vec4_D);CHKERRQ(ierr); 65 ierr = VecScatterDestroy(&pcbddc->R_to_B);CHKERRQ(ierr); 66 ierr = VecScatterDestroy(&pcbddc->R_to_D);CHKERRQ(ierr); 67 ierr = VecScatterDestroy(&pcbddc->coarse_loc_to_glob);CHKERRQ(ierr); 68 ierr = PetscFree(pcbddc->local_primal_indices);CHKERRQ(ierr); 69 ierr = PetscFree(pcbddc->replicated_local_primal_indices);CHKERRQ(ierr); 70 ierr = PetscFree(pcbddc->replicated_local_primal_values);CHKERRQ(ierr); 71 ierr = PetscFree(pcbddc->local_primal_displacements);CHKERRQ(ierr); 72 ierr = PetscFree(pcbddc->local_primal_sizes);CHKERRQ(ierr); 73 PetscFunctionReturn(0); 74 } 75 76 #undef __FUNCT__ 77 #define __FUNCT__ "PCBDDCCreateWorkVectors" 78 PetscErrorCode PCBDDCCreateWorkVectors(PC pc) 79 { 80 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 81 PC_IS *pcis = (PC_IS*)pc->data; 82 VecType impVecType; 83 PetscInt n_vertices,n_constraints,local_primal_size,n_R; 84 PetscErrorCode ierr; 85 86 PetscFunctionBegin; 87 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,NULL);CHKERRQ(ierr); 88 ierr = PCBDDCGetPrimalConstraintsLocalIdx(pc,&n_constraints,NULL);CHKERRQ(ierr); 89 local_primal_size = n_constraints+n_vertices; 90 n_R = pcis->n-n_vertices; 91 /* parallel work vectors used in presolve. TODO: move outside */ 92 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->original_rhs);CHKERRQ(ierr); 93 ierr = VecDuplicate(pcis->vec1_global,&pcbddc->temp_solution);CHKERRQ(ierr); 94 /* local work vectors */ 95 ierr = VecGetType(pcis->vec1_N,&impVecType);CHKERRQ(ierr); 96 ierr = VecDuplicate(pcis->vec1_D,&pcbddc->vec4_D);CHKERRQ(ierr); 97 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_R);CHKERRQ(ierr); 98 ierr = VecSetSizes(pcbddc->vec1_R,PETSC_DECIDE,n_R);CHKERRQ(ierr); 99 ierr = VecSetType(pcbddc->vec1_R,impVecType);CHKERRQ(ierr); 100 ierr = VecDuplicate(pcbddc->vec1_R,&pcbddc->vec2_R);CHKERRQ(ierr); 101 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_P);CHKERRQ(ierr); 102 ierr = VecSetSizes(pcbddc->vec1_P,PETSC_DECIDE,local_primal_size);CHKERRQ(ierr); 103 ierr = VecSetType(pcbddc->vec1_P,impVecType);CHKERRQ(ierr); 104 if (n_constraints) { 105 ierr = VecCreate(PetscObjectComm((PetscObject)pcis->vec1_N),&pcbddc->vec1_C);CHKERRQ(ierr); 106 ierr = VecSetSizes(pcbddc->vec1_C,PETSC_DECIDE,n_constraints);CHKERRQ(ierr); 107 ierr = VecSetType(pcbddc->vec1_C,impVecType);CHKERRQ(ierr); 108 } 109 PetscFunctionReturn(0); 110 } 111 112 #undef __FUNCT__ 113 #define __FUNCT__ "PCBDDCSetUpLocalMatrices" 114 PetscErrorCode PCBDDCSetUpLocalMatrices(PC pc) 115 { 116 PC_IS* pcis = (PC_IS*)(pc->data); 117 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 118 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 119 /* manage repeated solves */ 120 MatReuse reuse; 121 MatStructure matstruct; 122 PetscErrorCode ierr; 123 124 PetscFunctionBegin; 125 /* get mat flags */ 126 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 127 reuse = MAT_INITIAL_MATRIX; 128 if (pc->setupcalled) { 129 /* when matstruct is SAME_PRECONDITIONER, we shouldn't be here */ 130 if (matstruct == SAME_PRECONDITIONER) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"This should not happen"); 131 if (matstruct == SAME_NONZERO_PATTERN) { 132 reuse = MAT_REUSE_MATRIX; 133 } else { 134 reuse = MAT_INITIAL_MATRIX; 135 } 136 } 137 if (reuse == MAT_INITIAL_MATRIX) { 138 ierr = MatDestroy(&pcis->A_II);CHKERRQ(ierr); 139 ierr = MatDestroy(&pcis->A_IB);CHKERRQ(ierr); 140 ierr = MatDestroy(&pcis->A_BI);CHKERRQ(ierr); 141 ierr = MatDestroy(&pcis->A_BB);CHKERRQ(ierr); 142 ierr = MatDestroy(&pcbddc->local_mat);CHKERRQ(ierr); 143 } 144 145 /* transform local matrices if needed */ 146 if (pcbddc->use_change_of_basis) { 147 Mat change_mat_all; 148 PetscScalar *row_cmat_values; 149 PetscInt *row_cmat_indices; 150 PetscInt *nnz,*is_indices,*temp_indices; 151 PetscInt i,j,k,n_D,n_B; 152 153 /* Get Non-overlapping dimensions */ 154 n_B = pcis->n_B; 155 n_D = pcis->n-n_B; 156 157 /* compute nonzero structure of change of basis on all local nodes */ 158 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 159 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 160 for (i=0;i<n_D;i++) nnz[is_indices[i]] = 1; 161 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 162 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 163 k=1; 164 for (i=0;i<n_B;i++) { 165 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 166 nnz[is_indices[i]]=j; 167 if (k < j) k = j; 168 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,NULL,NULL);CHKERRQ(ierr); 169 } 170 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 171 /* assemble change of basis matrix on the whole set of local dofs */ 172 ierr = PetscMalloc(k*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 173 ierr = MatCreate(PETSC_COMM_SELF,&change_mat_all);CHKERRQ(ierr); 174 ierr = MatSetSizes(change_mat_all,pcis->n,pcis->n,pcis->n,pcis->n);CHKERRQ(ierr); 175 ierr = MatSetType(change_mat_all,MATSEQAIJ);CHKERRQ(ierr); 176 ierr = MatSeqAIJSetPreallocation(change_mat_all,0,nnz);CHKERRQ(ierr); 177 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 178 for (i=0;i<n_D;i++) { 179 ierr = MatSetValue(change_mat_all,is_indices[i],is_indices[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 180 } 181 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 182 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 183 for (i=0;i<n_B;i++) { 184 ierr = MatGetRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 185 for (k=0; k<j; k++) temp_indices[k]=is_indices[row_cmat_indices[k]]; 186 ierr = MatSetValues(change_mat_all,1,&is_indices[i],j,temp_indices,row_cmat_values,INSERT_VALUES);CHKERRQ(ierr); 187 ierr = MatRestoreRow(pcbddc->ChangeOfBasisMatrix,i,&j,(const PetscInt**)&row_cmat_indices,(const PetscScalar**)&row_cmat_values);CHKERRQ(ierr); 188 } 189 ierr = MatAssemblyBegin(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 190 ierr = MatAssemblyEnd(change_mat_all,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 191 /* TODO: HOW TO WORK WITH BAIJ? PtAP not provided */ 192 ierr = MatGetBlockSize(matis->A,&i);CHKERRQ(ierr); 193 if (i==1) { 194 ierr = MatPtAP(matis->A,change_mat_all,reuse,2.0,&pcbddc->local_mat);CHKERRQ(ierr); 195 } else { 196 Mat work_mat; 197 ierr = MatConvert(matis->A,MATSEQAIJ,MAT_INITIAL_MATRIX,&work_mat);CHKERRQ(ierr); 198 ierr = MatPtAP(work_mat,change_mat_all,reuse,2.0,&pcbddc->local_mat);CHKERRQ(ierr); 199 ierr = MatDestroy(&work_mat);CHKERRQ(ierr); 200 } 201 ierr = MatDestroy(&change_mat_all);CHKERRQ(ierr); 202 ierr = PetscFree(nnz);CHKERRQ(ierr); 203 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 204 } else { 205 /* without change of basis, the local matrix is unchanged */ 206 if (!pcbddc->local_mat) { 207 ierr = PetscObjectReference((PetscObject)matis->A);CHKERRQ(ierr); 208 pcbddc->local_mat = matis->A; 209 } 210 } 211 212 /* get submatrices */ 213 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_I_local,reuse,&pcis->A_II);CHKERRQ(ierr); 214 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_I_local,pcis->is_B_local,reuse,&pcis->A_IB);CHKERRQ(ierr); 215 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_I_local,reuse,&pcis->A_BI);CHKERRQ(ierr); 216 ierr = MatGetSubMatrix(pcbddc->local_mat,pcis->is_B_local,pcis->is_B_local,reuse,&pcis->A_BB);CHKERRQ(ierr); 217 PetscFunctionReturn(0); 218 } 219 220 #undef __FUNCT__ 221 #define __FUNCT__ "PCBDDCSetUpLocalScatters" 222 PetscErrorCode PCBDDCSetUpLocalScatters(PC pc,IS* is_R_local_n) 223 { 224 PC_IS* pcis = (PC_IS*)(pc->data); 225 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 226 IS is_R_local,is_aux1,is_aux2; 227 PetscInt *vertices,*aux_array1,*aux_array2,*is_indices,*idx_R_local; 228 PetscInt n_vertices,n_constraints,i,j,n_R,n_D,n_B; 229 PetscBool *array_bool; 230 PetscErrorCode ierr; 231 232 PetscFunctionBegin; 233 /* Set Non-overlapping dimensions */ 234 n_B = pcis->n_B; n_D = pcis->n - n_B; 235 /* get vertex indices from constraint matrix */ 236 ierr = PCBDDCGetPrimalVerticesLocalIdx(pc,&n_vertices,&vertices);CHKERRQ(ierr); 237 /* Set number of constraints */ 238 n_constraints = pcbddc->local_primal_size-n_vertices; 239 /* Dohrmann's notation: dofs splitted in R (Remaining: all dofs but the vertices) and V (Vertices) */ 240 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&array_bool);CHKERRQ(ierr); 241 for (i=0;i<pcis->n;i++) array_bool[i] = PETSC_TRUE; 242 for (i=0;i<n_vertices;i++) array_bool[vertices[i]] = PETSC_FALSE; 243 ierr = PetscMalloc((pcis->n-n_vertices)*sizeof(PetscInt),&idx_R_local);CHKERRQ(ierr); 244 for (i=0, n_R=0; i<pcis->n; i++) { 245 if (array_bool[i]) { 246 idx_R_local[n_R] = i; 247 n_R++; 248 } 249 } 250 ierr = PetscFree(vertices);CHKERRQ(ierr); 251 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_R,idx_R_local,PETSC_OWN_POINTER,&is_R_local);CHKERRQ(ierr); 252 253 /* print some info if requested */ 254 if (pcbddc->dbg_flag) { 255 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 256 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 257 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d local dimensions\n",PetscGlobalRank);CHKERRQ(ierr); 258 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"local_size = %d, dirichlet_size = %d, boundary_size = %d\n",pcis->n,n_D,n_B);CHKERRQ(ierr); 259 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); 260 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"pcbddc->n_vertices = %d, pcbddc->n_constraints = %d\n",pcbddc->n_vertices,pcbddc->n_constraints);CHKERRQ(ierr); 261 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 262 } 263 264 /* VecScatters pcbddc->R_to_B and (optionally) pcbddc->R_to_D */ 265 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 266 ierr = PetscMalloc((pcis->n_B-n_vertices)*sizeof(PetscInt),&aux_array2);CHKERRQ(ierr); 267 ierr = ISGetIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 268 for (i=0; i<n_D; i++) array_bool[is_indices[i]] = PETSC_FALSE; 269 ierr = ISRestoreIndices(pcis->is_I_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 270 for (i=0, j=0; i<n_R; i++) { 271 if (array_bool[idx_R_local[i]]) { 272 aux_array1[j] = i; 273 j++; 274 } 275 } 276 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_OWN_POINTER,&is_aux1);CHKERRQ(ierr); 277 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 278 for (i=0, j=0; i<n_B; i++) { 279 if (array_bool[is_indices[i]]) { 280 aux_array2[j] = i; j++; 281 } 282 } 283 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 284 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array2,PETSC_OWN_POINTER,&is_aux2);CHKERRQ(ierr); 285 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_B,is_aux2,&pcbddc->R_to_B);CHKERRQ(ierr); 286 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 287 ierr = ISDestroy(&is_aux2);CHKERRQ(ierr); 288 289 if (pcbddc->inexact_prec_type || pcbddc->dbg_flag ) { 290 ierr = PetscMalloc(n_D*sizeof(PetscInt),&aux_array1);CHKERRQ(ierr); 291 for (i=0, j=0; i<n_R; i++) { 292 if (!array_bool[idx_R_local[i]]) { 293 aux_array1[j] = i; 294 j++; 295 } 296 } 297 ierr = ISCreateGeneral(PETSC_COMM_SELF,j,aux_array1,PETSC_OWN_POINTER,&is_aux1);CHKERRQ(ierr); 298 ierr = VecScatterCreate(pcbddc->vec1_R,is_aux1,pcis->vec1_D,(IS)0,&pcbddc->R_to_D);CHKERRQ(ierr); 299 ierr = ISDestroy(&is_aux1);CHKERRQ(ierr); 300 } 301 ierr = PetscFree(array_bool);CHKERRQ(ierr); 302 *is_R_local_n = is_R_local; 303 PetscFunctionReturn(0); 304 } 305 306 #undef __FUNCT__ 307 #define __FUNCT__ "PCBDDCSetUseExactDirichlet" 308 PetscErrorCode PCBDDCSetUseExactDirichlet(PC pc,PetscBool use) 309 { 310 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 311 312 PetscFunctionBegin; 313 pcbddc->use_exact_dirichlet=use; 314 PetscFunctionReturn(0); 315 } 316 317 #undef __FUNCT__ 318 #define __FUNCT__ "PCBDDCSetUpLocalSolvers" 319 PetscErrorCode PCBDDCSetUpLocalSolvers(PC pc, IS is_I_local, IS is_R_local) 320 { 321 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 322 PC_IS *pcis = (PC_IS*)pc->data; 323 PC pc_temp; 324 Mat A_RR; 325 Vec vec1,vec2,vec3; 326 MatStructure matstruct; 327 PetscScalar m_one = -1.0; 328 PetscReal value; 329 PetscInt n_D,n_R,use_exact,use_exact_reduced; 330 PetscErrorCode ierr; 331 332 PetscFunctionBegin; 333 /* Creating PC contexts for local Dirichlet and Neumann problems */ 334 ierr = PCGetOperators(pc,NULL,NULL,&matstruct);CHKERRQ(ierr); 335 336 /* DIRICHLET PROBLEM */ 337 /* Matrix for Dirichlet problem is pcis->A_II */ 338 ierr = ISGetSize(is_I_local,&n_D);CHKERRQ(ierr); 339 if (!pcbddc->ksp_D) { /* create object if not yet build */ 340 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_D);CHKERRQ(ierr); 341 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_D,(PetscObject)pc,1);CHKERRQ(ierr); 342 /* default */ 343 ierr = KSPSetType(pcbddc->ksp_D,KSPPREONLY);CHKERRQ(ierr); 344 ierr = KSPSetOptionsPrefix(pcbddc->ksp_D,"dirichlet_");CHKERRQ(ierr); 345 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 346 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 347 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 348 } 349 ierr = KSPSetOperators(pcbddc->ksp_D,pcis->A_II,pcis->A_II,matstruct);CHKERRQ(ierr); 350 /* Allow user's customization */ 351 ierr = KSPSetFromOptions(pcbddc->ksp_D);CHKERRQ(ierr); 352 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 353 if (!n_D) { 354 ierr = KSPGetPC(pcbddc->ksp_D,&pc_temp);CHKERRQ(ierr); 355 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 356 } 357 /* Set Up KSP for Dirichlet problem of BDDC */ 358 ierr = KSPSetUp(pcbddc->ksp_D);CHKERRQ(ierr); 359 /* set ksp_D into pcis data */ 360 ierr = KSPDestroy(&pcis->ksp_D);CHKERRQ(ierr); 361 ierr = PetscObjectReference((PetscObject)pcbddc->ksp_D);CHKERRQ(ierr); 362 pcis->ksp_D = pcbddc->ksp_D; 363 364 /* NEUMANN PROBLEM */ 365 /* Matrix for Neumann problem is A_RR -> we need to create it */ 366 ierr = ISGetSize(is_R_local,&n_R);CHKERRQ(ierr); 367 ierr = MatGetSubMatrix(pcbddc->local_mat,is_R_local,is_R_local,MAT_INITIAL_MATRIX,&A_RR);CHKERRQ(ierr); 368 if (!pcbddc->ksp_R) { /* create object if not yet build */ 369 ierr = KSPCreate(PETSC_COMM_SELF,&pcbddc->ksp_R);CHKERRQ(ierr); 370 ierr = PetscObjectIncrementTabLevel((PetscObject)pcbddc->ksp_R,(PetscObject)pc,1);CHKERRQ(ierr); 371 /* default */ 372 ierr = KSPSetType(pcbddc->ksp_R,KSPPREONLY);CHKERRQ(ierr); 373 ierr = KSPSetOptionsPrefix(pcbddc->ksp_R,"neumann_");CHKERRQ(ierr); 374 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 375 ierr = PCSetType(pc_temp,PCLU);CHKERRQ(ierr); 376 ierr = PCFactorSetReuseFill(pc_temp,PETSC_TRUE);CHKERRQ(ierr); 377 } 378 ierr = KSPSetOperators(pcbddc->ksp_R,A_RR,A_RR,matstruct);CHKERRQ(ierr); 379 /* Allow user's customization */ 380 ierr = KSPSetFromOptions(pcbddc->ksp_R);CHKERRQ(ierr); 381 /* umfpack interface has a bug when matrix dimension is zero. TODO solve from umfpack interface */ 382 if (!n_R) { 383 ierr = KSPGetPC(pcbddc->ksp_R,&pc_temp);CHKERRQ(ierr); 384 ierr = PCSetType(pc_temp,PCNONE);CHKERRQ(ierr); 385 } 386 /* Set Up KSP for Neumann problem of BDDC */ 387 ierr = KSPSetUp(pcbddc->ksp_R);CHKERRQ(ierr); 388 389 /* check Dirichlet and Neumann solvers and adapt them if a nullspace correction is needed */ 390 391 /* Dirichlet */ 392 ierr = MatGetVecs(pcis->A_II,&vec1,&vec2);CHKERRQ(ierr); 393 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 394 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 395 ierr = MatMult(pcis->A_II,vec1,vec2);CHKERRQ(ierr); 396 ierr = KSPSolve(pcbddc->ksp_D,vec2,vec3);CHKERRQ(ierr); 397 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 398 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 399 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 400 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 401 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 402 /* need to be adapted? */ 403 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 404 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 405 ierr = PCBDDCSetUseExactDirichlet(pc,(PetscBool)use_exact_reduced);CHKERRQ(ierr); 406 /* print info */ 407 if (pcbddc->dbg_flag) { 408 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 409 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"--------------------------------------------------\n");CHKERRQ(ierr); 410 ierr = PetscViewerASCIIPrintf(pcbddc->dbg_viewer,"Checking solution of Dirichlet and Neumann problems\n");CHKERRQ(ierr); 411 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Dirichlet solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 412 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 413 } 414 if (n_D && pcbddc->NullSpace && !use_exact_reduced && !pcbddc->inexact_prec_type) { 415 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_I_local);CHKERRQ(ierr); 416 } 417 418 /* Neumann */ 419 ierr = MatGetVecs(A_RR,&vec1,&vec2);CHKERRQ(ierr); 420 ierr = VecDuplicate(vec1,&vec3);CHKERRQ(ierr); 421 ierr = VecSetRandom(vec1,NULL);CHKERRQ(ierr); 422 ierr = MatMult(A_RR,vec1,vec2);CHKERRQ(ierr); 423 ierr = KSPSolve(pcbddc->ksp_R,vec2,vec3);CHKERRQ(ierr); 424 ierr = VecAXPY(vec3,m_one,vec1);CHKERRQ(ierr); 425 ierr = VecNorm(vec3,NORM_INFINITY,&value);CHKERRQ(ierr); 426 ierr = VecDestroy(&vec1);CHKERRQ(ierr); 427 ierr = VecDestroy(&vec2);CHKERRQ(ierr); 428 ierr = VecDestroy(&vec3);CHKERRQ(ierr); 429 /* need to be adapted? */ 430 use_exact = (PetscAbsReal(value) > 1.e-4 ? 0 : 1); 431 if (PetscAbsReal(value) > 1.e-4) use_exact = 0; 432 ierr = MPI_Allreduce(&use_exact,&use_exact_reduced,1,MPIU_INT,MPI_LAND,PetscObjectComm((PetscObject)pc));CHKERRQ(ierr); 433 /* print info */ 434 if (pcbddc->dbg_flag) { 435 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d infinity error for Neumann solve = % 1.14e \n",PetscGlobalRank,value);CHKERRQ(ierr); 436 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 437 } 438 if (n_R && pcbddc->NullSpace && !use_exact_reduced) { /* is it the right logic? */ 439 ierr = PCBDDCNullSpaceAssembleCorrection(pc,is_R_local);CHKERRQ(ierr); 440 } 441 442 /* free Neumann problem's matrix */ 443 ierr = MatDestroy(&A_RR);CHKERRQ(ierr); 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNCT__ 448 #define __FUNCT__ "PCBDDCSolveSaddlePoint" 449 static PetscErrorCode PCBDDCSolveSaddlePoint(PC pc) 450 { 451 PetscErrorCode ierr; 452 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 453 454 PetscFunctionBegin; 455 ierr = KSPSolve(pcbddc->ksp_R,pcbddc->vec1_R,pcbddc->vec2_R);CHKERRQ(ierr); 456 if (pcbddc->local_auxmat1) { 457 ierr = MatMult(pcbddc->local_auxmat1,pcbddc->vec2_R,pcbddc->vec1_C);CHKERRQ(ierr); 458 ierr = MatMultAdd(pcbddc->local_auxmat2,pcbddc->vec1_C,pcbddc->vec2_R,pcbddc->vec2_R);CHKERRQ(ierr); 459 } 460 PetscFunctionReturn(0); 461 } 462 463 #undef __FUNCT__ 464 #define __FUNCT__ "PCBDDCApplyInterfacePreconditioner" 465 PetscErrorCode PCBDDCApplyInterfacePreconditioner(PC pc) 466 { 467 PetscErrorCode ierr; 468 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 469 PC_IS* pcis = (PC_IS*) (pc->data); 470 const PetscScalar zero = 0.0; 471 472 PetscFunctionBegin; 473 /* Application of PHI^T (or PSI^T) */ 474 if (pcbddc->coarse_psi_B) { 475 ierr = MatMultTranspose(pcbddc->coarse_psi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 476 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_psi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 477 } else { 478 ierr = MatMultTranspose(pcbddc->coarse_phi_B,pcis->vec1_B,pcbddc->vec1_P);CHKERRQ(ierr); 479 if (pcbddc->inexact_prec_type) { ierr = MatMultTransposeAdd(pcbddc->coarse_phi_D,pcis->vec1_D,pcbddc->vec1_P,pcbddc->vec1_P);CHKERRQ(ierr); } 480 } 481 /* Scatter data of coarse_rhs */ 482 if (pcbddc->coarse_rhs) { ierr = VecSet(pcbddc->coarse_rhs,zero);CHKERRQ(ierr); } 483 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 484 485 /* Local solution on R nodes */ 486 ierr = VecSet(pcbddc->vec1_R,zero);CHKERRQ(ierr); 487 ierr = VecScatterBegin(pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 488 ierr = VecScatterEnd (pcbddc->R_to_B,pcis->vec1_B,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 489 if (pcbddc->inexact_prec_type) { 490 ierr = VecScatterBegin(pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 491 ierr = VecScatterEnd (pcbddc->R_to_D,pcis->vec1_D,pcbddc->vec1_R,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 492 } 493 ierr = PCBDDCSolveSaddlePoint(pc);CHKERRQ(ierr); 494 ierr = VecSet(pcis->vec1_B,zero);CHKERRQ(ierr); 495 ierr = VecScatterBegin(pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 496 ierr = VecScatterEnd (pcbddc->R_to_B,pcbddc->vec2_R,pcis->vec1_B,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 497 if (pcbddc->inexact_prec_type) { 498 ierr = VecScatterBegin(pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 499 ierr = VecScatterEnd (pcbddc->R_to_D,pcbddc->vec2_R,pcis->vec1_D,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 500 } 501 502 /* Coarse solution */ 503 ierr = PCBDDCScatterCoarseDataEnd(pc,pcbddc->vec1_P,pcbddc->coarse_rhs,ADD_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 504 if (pcbddc->coarse_rhs) { /* TODO remove null space when doing multilevel */ 505 ierr = KSPSolve(pcbddc->coarse_ksp,pcbddc->coarse_rhs,pcbddc->coarse_vec);CHKERRQ(ierr); 506 } 507 ierr = PCBDDCScatterCoarseDataBegin(pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 508 ierr = PCBDDCScatterCoarseDataEnd (pc,pcbddc->coarse_vec,pcbddc->vec1_P,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 509 510 /* Sum contributions from two levels */ 511 ierr = MatMultAdd(pcbddc->coarse_phi_B,pcbddc->vec1_P,pcis->vec1_B,pcis->vec1_B);CHKERRQ(ierr); 512 if (pcbddc->inexact_prec_type) { ierr = MatMultAdd(pcbddc->coarse_phi_D,pcbddc->vec1_P,pcis->vec1_D,pcis->vec1_D);CHKERRQ(ierr); } 513 PetscFunctionReturn(0); 514 } 515 516 #undef __FUNCT__ 517 #define __FUNCT__ "PCBDDCScatterCoarseDataBegin" 518 PetscErrorCode PCBDDCScatterCoarseDataBegin(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 519 { 520 PetscErrorCode ierr; 521 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 522 523 PetscFunctionBegin; 524 switch (pcbddc->coarse_communications_type) { 525 case SCATTERS_BDDC: 526 ierr = VecScatterBegin(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 527 break; 528 case GATHERS_BDDC: 529 break; 530 } 531 PetscFunctionReturn(0); 532 } 533 534 #undef __FUNCT__ 535 #define __FUNCT__ "PCBDDCScatterCoarseDataEnd" 536 PetscErrorCode PCBDDCScatterCoarseDataEnd(PC pc,Vec vec_from, Vec vec_to, InsertMode imode, ScatterMode smode) 537 { 538 PetscErrorCode ierr; 539 PC_BDDC* pcbddc = (PC_BDDC*)(pc->data); 540 PetscScalar* array_to; 541 PetscScalar* array_from; 542 MPI_Comm comm; 543 PetscInt i; 544 545 PetscFunctionBegin; 546 ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr); 547 switch (pcbddc->coarse_communications_type) { 548 case SCATTERS_BDDC: 549 ierr = VecScatterEnd(pcbddc->coarse_loc_to_glob,vec_from,vec_to,imode,smode);CHKERRQ(ierr); 550 break; 551 case GATHERS_BDDC: 552 if (vec_from) { 553 ierr = VecGetArray(vec_from,&array_from);CHKERRQ(ierr); 554 } 555 if (vec_to) { 556 ierr = VecGetArray(vec_to,&array_to);CHKERRQ(ierr); 557 } 558 switch(pcbddc->coarse_problem_type){ 559 case SEQUENTIAL_BDDC: 560 if (smode == SCATTER_FORWARD) { 561 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); 562 if (vec_to) { 563 if (imode == ADD_VALUES) { 564 for (i=0;i<pcbddc->replicated_primal_size;i++) { 565 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 566 } 567 } else { 568 for (i=0;i<pcbddc->replicated_primal_size;i++) { 569 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 570 } 571 } 572 } 573 } else { 574 if (vec_from) { 575 if (imode == ADD_VALUES) { 576 MPI_Comm vec_from_comm; 577 ierr = PetscObjectGetComm((PetscObject)(vec_from),&vec_from_comm);CHKERRQ(ierr); 578 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); 579 } 580 for (i=0;i<pcbddc->replicated_primal_size;i++) { 581 pcbddc->replicated_local_primal_values[i]=array_from[pcbddc->replicated_local_primal_indices[i]]; 582 } 583 } 584 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); 585 } 586 break; 587 case REPLICATED_BDDC: 588 if (smode == SCATTER_FORWARD) { 589 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); 590 if (imode == ADD_VALUES) { 591 for (i=0;i<pcbddc->replicated_primal_size;i++) { 592 array_to[pcbddc->replicated_local_primal_indices[i]]+=pcbddc->replicated_local_primal_values[i]; 593 } 594 } else { 595 for (i=0;i<pcbddc->replicated_primal_size;i++) { 596 array_to[pcbddc->replicated_local_primal_indices[i]]=pcbddc->replicated_local_primal_values[i]; 597 } 598 } 599 } else { /* no communications needed for SCATTER_REVERSE since needed data is already present */ 600 if (imode == ADD_VALUES) { 601 for (i=0;i<pcbddc->local_primal_size;i++) { 602 array_to[i]+=array_from[pcbddc->local_primal_indices[i]]; 603 } 604 } else { 605 for (i=0;i<pcbddc->local_primal_size;i++) { 606 array_to[i]=array_from[pcbddc->local_primal_indices[i]]; 607 } 608 } 609 } 610 break; 611 case MULTILEVEL_BDDC: 612 break; 613 case PARALLEL_BDDC: 614 break; 615 } 616 if (vec_from) { 617 ierr = VecRestoreArray(vec_from,&array_from);CHKERRQ(ierr); 618 } 619 if (vec_to) { 620 ierr = VecRestoreArray(vec_to,&array_to);CHKERRQ(ierr); 621 } 622 break; 623 } 624 PetscFunctionReturn(0); 625 } 626 627 /* uncomment for testing purposes */ 628 /* #define PETSC_MISSING_LAPACK_GESVD 1 */ 629 #undef __FUNCT__ 630 #define __FUNCT__ "PCBDDCConstraintsSetUp" 631 PetscErrorCode PCBDDCConstraintsSetUp(PC pc) 632 { 633 PetscErrorCode ierr; 634 PC_IS* pcis = (PC_IS*)(pc->data); 635 PC_BDDC* pcbddc = (PC_BDDC*)pc->data; 636 Mat_IS* matis = (Mat_IS*)pc->pmat->data; 637 /* constraint and (optionally) change of basis matrix implemented as SeqAIJ */ 638 MatType impMatType=MATSEQAIJ; 639 /* one and zero */ 640 PetscScalar one=1.0,zero=0.0; 641 /* space to store constraints and their local indices */ 642 PetscScalar *temp_quadrature_constraint; 643 PetscInt *temp_indices,*temp_indices_to_constraint,*temp_indices_to_constraint_B; 644 /* iterators */ 645 PetscInt i,j,k,total_counts,temp_start_ptr; 646 /* stuff to store connected components stored in pcbddc->mat_graph */ 647 IS ISForVertices,*ISForFaces,*ISForEdges,*used_IS; 648 PetscInt n_ISForFaces,n_ISForEdges; 649 PetscBool get_faces,get_edges,get_vertices; 650 /* near null space stuff */ 651 MatNullSpace nearnullsp; 652 const Vec *nearnullvecs; 653 Vec *localnearnullsp; 654 PetscBool nnsp_has_cnst; 655 PetscInt nnsp_size; 656 PetscScalar *array; 657 /* BLAS integers */ 658 PetscBLASInt lwork,lierr; 659 PetscBLASInt Blas_N,Blas_M,Blas_K,Blas_one=1; 660 PetscBLASInt Blas_LDA,Blas_LDB,Blas_LDC; 661 /* LAPACK working arrays for SVD or POD */ 662 PetscBool skip_lapack; 663 PetscScalar *work; 664 PetscReal *singular_vals; 665 #if defined(PETSC_USE_COMPLEX) 666 PetscReal *rwork; 667 #endif 668 #if defined(PETSC_MISSING_LAPACK_GESVD) 669 PetscBLASInt Blas_one_2=1; 670 PetscScalar *temp_basis,*correlation_mat; 671 #else 672 PetscBLASInt dummy_int_1=1,dummy_int_2=1; 673 PetscScalar dummy_scalar_1=0.0,dummy_scalar_2=0.0; 674 #endif 675 /* change of basis */ 676 PetscInt *aux_primal_numbering,*aux_primal_minloc,*global_indices; 677 PetscBool boolforchange,*change_basis,*touched; 678 /* auxiliary stuff */ 679 PetscInt *nnz,*is_indices,*local_to_B; 680 /* some quantities */ 681 PetscInt n_vertices,total_primal_vertices; 682 PetscInt size_of_constraint,max_size_of_constraint,max_constraints,temp_constraints; 683 684 685 PetscFunctionBegin; 686 /* Get index sets for faces, edges and vertices from graph */ 687 get_faces = PETSC_TRUE; 688 get_edges = PETSC_TRUE; 689 get_vertices = PETSC_TRUE; 690 if (pcbddc->vertices_flag) { 691 get_faces = PETSC_FALSE; 692 get_edges = PETSC_FALSE; 693 } 694 if (pcbddc->constraints_flag) { 695 get_vertices = PETSC_FALSE; 696 } 697 if (pcbddc->faces_flag) { 698 get_edges = PETSC_FALSE; 699 } 700 if (pcbddc->edges_flag) { 701 get_faces = PETSC_FALSE; 702 } 703 /* default */ 704 if (!get_faces && !get_edges && !get_vertices) { 705 get_vertices = PETSC_TRUE; 706 } 707 ierr = PCBDDCGraphGetCandidatesIS(pcbddc->mat_graph,get_faces,get_edges,get_vertices,&n_ISForFaces,&ISForFaces,&n_ISForEdges,&ISForEdges,&ISForVertices); 708 /* print some info */ 709 if (pcbddc->dbg_flag) { 710 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 711 i = 0; 712 if (ISForVertices) { 713 ierr = ISGetSize(ISForVertices,&i);CHKERRQ(ierr); 714 } 715 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate vertices\n",PetscGlobalRank,i);CHKERRQ(ierr); 716 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate edges\n",PetscGlobalRank,n_ISForEdges);CHKERRQ(ierr); 717 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Subdomain %04d got %02d local candidate faces\n",PetscGlobalRank,n_ISForFaces);CHKERRQ(ierr); 718 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 719 } 720 /* check if near null space is attached to global mat */ 721 ierr = MatGetNearNullSpace(pc->pmat,&nearnullsp);CHKERRQ(ierr); 722 if (nearnullsp) { 723 ierr = MatNullSpaceGetVecs(nearnullsp,&nnsp_has_cnst,&nnsp_size,&nearnullvecs);CHKERRQ(ierr); 724 } else { /* if near null space is not provided BDDC uses constants by default */ 725 nnsp_size = 0; 726 nnsp_has_cnst = PETSC_TRUE; 727 } 728 /* get max number of constraints on a single cc */ 729 max_constraints = nnsp_size; 730 if (nnsp_has_cnst) max_constraints++; 731 732 /* 733 Evaluate maximum storage size needed by the procedure 734 - temp_indices will contain start index of each constraint stored as follows 735 - temp_indices_to_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the indices (in local numbering) on which the constraint acts 736 - 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 737 - temp_quadrature_constraint [temp_indices[i],...,temp[indices[i+1]-1] will contain the scalars representing the constraint itself 738 */ 739 total_counts = n_ISForFaces+n_ISForEdges; 740 total_counts *= max_constraints; 741 n_vertices = 0; 742 if (ISForVertices) { 743 ierr = ISGetSize(ISForVertices,&n_vertices);CHKERRQ(ierr); 744 } 745 total_counts += n_vertices; 746 ierr = PetscMalloc((total_counts+1)*sizeof(PetscInt),&temp_indices);CHKERRQ(ierr); 747 ierr = PetscMalloc((total_counts+1)*sizeof(PetscBool),&change_basis);CHKERRQ(ierr); 748 total_counts = 0; 749 max_size_of_constraint = 0; 750 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 751 if (i<n_ISForEdges) { 752 used_IS = &ISForEdges[i]; 753 } else { 754 used_IS = &ISForFaces[i-n_ISForEdges]; 755 } 756 ierr = ISGetSize(*used_IS,&j);CHKERRQ(ierr); 757 total_counts += j; 758 max_size_of_constraint = PetscMax(j,max_size_of_constraint); 759 } 760 total_counts *= max_constraints; 761 total_counts += n_vertices; 762 ierr = PetscMalloc(total_counts*sizeof(PetscScalar),&temp_quadrature_constraint);CHKERRQ(ierr); 763 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint);CHKERRQ(ierr); 764 ierr = PetscMalloc(total_counts*sizeof(PetscInt),&temp_indices_to_constraint_B);CHKERRQ(ierr); 765 /* local to boundary numbering */ 766 ierr = PetscMalloc(pcis->n*sizeof(PetscInt),&local_to_B);CHKERRQ(ierr); 767 ierr = ISGetIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 768 for (i=0;i<pcis->n;i++) local_to_B[i]=-1; 769 for (i=0;i<pcis->n_B;i++) local_to_B[is_indices[i]]=i; 770 ierr = ISRestoreIndices(pcis->is_B_local,(const PetscInt**)&is_indices);CHKERRQ(ierr); 771 /* get local part of global near null space vectors */ 772 ierr = PetscMalloc(nnsp_size*sizeof(Vec),&localnearnullsp);CHKERRQ(ierr); 773 for (k=0;k<nnsp_size;k++) { 774 ierr = VecDuplicate(pcis->vec1_N,&localnearnullsp[k]);CHKERRQ(ierr); 775 ierr = VecScatterBegin(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 776 ierr = VecScatterEnd(matis->ctx,nearnullvecs[k],localnearnullsp[k],INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 777 } 778 779 /* whether or not to skip lapack calls */ 780 skip_lapack = PETSC_TRUE; 781 if (n_ISForFaces+n_ISForEdges) skip_lapack = PETSC_FALSE; 782 783 /* First we issue queries to allocate optimal workspace for LAPACKgesvd (or LAPACKsyev if SVD is missing) */ 784 if (!pcbddc->use_nnsp_true && !skip_lapack) { 785 PetscScalar temp_work; 786 #if defined(PETSC_MISSING_LAPACK_GESVD) 787 /* Proper Orthogonal Decomposition (POD) using the snapshot method */ 788 ierr = PetscMalloc(max_constraints*max_constraints*sizeof(PetscScalar),&correlation_mat);CHKERRQ(ierr); 789 ierr = PetscMalloc(max_constraints*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 790 ierr = PetscMalloc(max_size_of_constraint*max_constraints*sizeof(PetscScalar),&temp_basis);CHKERRQ(ierr); 791 #if defined(PETSC_USE_COMPLEX) 792 ierr = PetscMalloc(3*max_constraints*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 793 #endif 794 /* now we evaluate the optimal workspace using query with lwork=-1 */ 795 ierr = PetscBLASIntCast(max_constraints,&Blas_N);CHKERRQ(ierr); 796 ierr = PetscBLASIntCast(max_constraints,&Blas_LDA);CHKERRQ(ierr); 797 lwork = -1; 798 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 799 #if !defined(PETSC_USE_COMPLEX) 800 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,&temp_work,&lwork,&lierr)); 801 #else 802 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,&temp_work,&lwork,rwork,&lierr)); 803 #endif 804 ierr = PetscFPTrapPop();CHKERRQ(ierr); 805 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to SYEV Lapack routine %d",(int)lierr); 806 #else /* on missing GESVD */ 807 /* SVD */ 808 PetscInt max_n,min_n; 809 max_n = max_size_of_constraint; 810 min_n = max_constraints; 811 if (max_size_of_constraint < max_constraints) { 812 min_n = max_size_of_constraint; 813 max_n = max_constraints; 814 } 815 ierr = PetscMalloc(min_n*sizeof(PetscReal),&singular_vals);CHKERRQ(ierr); 816 #if defined(PETSC_USE_COMPLEX) 817 ierr = PetscMalloc(5*min_n*sizeof(PetscReal),&rwork);CHKERRQ(ierr); 818 #endif 819 /* now we evaluate the optimal workspace using query with lwork=-1 */ 820 lwork = -1; 821 ierr = PetscBLASIntCast(max_n,&Blas_M);CHKERRQ(ierr); 822 ierr = PetscBLASIntCast(min_n,&Blas_N);CHKERRQ(ierr); 823 ierr = PetscBLASIntCast(max_n,&Blas_LDA);CHKERRQ(ierr); 824 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 825 #if !defined(PETSC_USE_COMPLEX) 826 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)); 827 #else 828 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)); 829 #endif 830 ierr = PetscFPTrapPop();CHKERRQ(ierr); 831 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GESVD Lapack routine %d",(int)lierr); 832 #endif /* on missing GESVD */ 833 /* Allocate optimal workspace */ 834 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(temp_work),&lwork);CHKERRQ(ierr); 835 ierr = PetscMalloc((PetscInt)lwork*sizeof(PetscScalar),&work);CHKERRQ(ierr); 836 } 837 /* Now we can loop on constraining sets */ 838 total_counts = 0; 839 temp_indices[0] = 0; 840 /* vertices */ 841 if (ISForVertices) { 842 ierr = ISGetIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 843 if (nnsp_has_cnst) { /* consider all vertices */ 844 for (i=0;i<n_vertices;i++) { 845 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 846 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 847 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 848 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 849 change_basis[total_counts]=PETSC_FALSE; 850 total_counts++; 851 } 852 } else { /* consider vertices for which exist at least a localnearnullsp which is not null there */ 853 PetscBool used_vertex; 854 for (i=0;i<n_vertices;i++) { 855 used_vertex = PETSC_FALSE; 856 k = 0; 857 while (!used_vertex && k<nnsp_size) { 858 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 859 if (PetscAbsScalar(array[is_indices[i]])>0.0) { 860 temp_indices_to_constraint[temp_indices[total_counts]]=is_indices[i]; 861 temp_indices_to_constraint_B[temp_indices[total_counts]]=local_to_B[is_indices[i]]; 862 temp_quadrature_constraint[temp_indices[total_counts]]=1.0; 863 temp_indices[total_counts+1]=temp_indices[total_counts]+1; 864 change_basis[total_counts]=PETSC_FALSE; 865 total_counts++; 866 used_vertex = PETSC_TRUE; 867 } 868 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 869 k++; 870 } 871 } 872 } 873 ierr = ISRestoreIndices(ISForVertices,(const PetscInt**)&is_indices);CHKERRQ(ierr); 874 n_vertices = total_counts; 875 } 876 877 /* edges and faces */ 878 for (i=0;i<n_ISForEdges+n_ISForFaces;i++) { 879 if (i<n_ISForEdges) { 880 used_IS = &ISForEdges[i]; 881 boolforchange = pcbddc->use_change_of_basis; /* change or not the basis on the edge */ 882 } else { 883 used_IS = &ISForFaces[i-n_ISForEdges]; 884 boolforchange = (PetscBool)(pcbddc->use_change_of_basis && pcbddc->use_change_on_faces); /* change or not the basis on the face */ 885 } 886 temp_constraints = 0; /* zero the number of constraints I have on this conn comp */ 887 temp_start_ptr = total_counts; /* need to know the starting index of constraints stored */ 888 ierr = ISGetSize(*used_IS,&size_of_constraint);CHKERRQ(ierr); 889 ierr = ISGetIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 890 /* change of basis should not be performed on local periodic nodes */ 891 if (pcbddc->mat_graph->mirrors && pcbddc->mat_graph->mirrors[is_indices[0]]) boolforchange = PETSC_FALSE; 892 if (nnsp_has_cnst) { 893 PetscScalar quad_value; 894 temp_constraints++; 895 quad_value = (PetscScalar)(1.0/PetscSqrtReal((PetscReal)size_of_constraint)); 896 for (j=0;j<size_of_constraint;j++) { 897 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 898 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 899 temp_quadrature_constraint[temp_indices[total_counts]+j]=quad_value; 900 } 901 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 902 change_basis[total_counts]=boolforchange; 903 total_counts++; 904 } 905 for (k=0;k<nnsp_size;k++) { 906 PetscReal real_value; 907 ierr = VecGetArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 908 for (j=0;j<size_of_constraint;j++) { 909 temp_indices_to_constraint[temp_indices[total_counts]+j]=is_indices[j]; 910 temp_indices_to_constraint_B[temp_indices[total_counts]+j]=local_to_B[is_indices[j]]; 911 temp_quadrature_constraint[temp_indices[total_counts]+j]=array[is_indices[j]]; 912 } 913 ierr = VecRestoreArrayRead(localnearnullsp[k],(const PetscScalar**)&array);CHKERRQ(ierr); 914 /* check if array is null on the connected component */ 915 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 916 PetscStackCallBLAS("BLASasum",real_value = BLASasum_(&Blas_N,&temp_quadrature_constraint[temp_indices[total_counts]],&Blas_one)); 917 if (real_value > 0.0) { /* keep indices and values */ 918 temp_constraints++; 919 temp_indices[total_counts+1]=temp_indices[total_counts]+size_of_constraint; /* store new starting point */ 920 change_basis[total_counts]=boolforchange; 921 total_counts++; 922 } 923 } 924 ierr = ISRestoreIndices(*used_IS,(const PetscInt**)&is_indices);CHKERRQ(ierr); 925 /* perform SVD on the constraints if use_nnsp_true has not be requested by the user */ 926 if (!pcbddc->use_nnsp_true) { 927 PetscReal tol = 1.0e-8; /* tolerance for retaining eigenmodes */ 928 929 #if defined(PETSC_MISSING_LAPACK_GESVD) 930 /* SVD: Y = U*S*V^H -> U (eigenvectors of Y*Y^H) = Y*V*(S)^\dag 931 POD: Y^H*Y = V*D*V^H, D = S^H*S -> U = Y*V*D^(-1/2) 932 -> When PETSC_USE_COMPLEX and PETSC_MISSING_LAPACK_GESVD are defined 933 the constraints basis will differ (by a complex factor with absolute value equal to 1) 934 from that computed using LAPACKgesvd 935 -> This is due to a different computation of eigenvectors in LAPACKheev 936 -> The quality of the POD-computed basis will be the same */ 937 ierr = PetscMemzero(correlation_mat,temp_constraints*temp_constraints*sizeof(PetscScalar));CHKERRQ(ierr); 938 /* Store upper triangular part of correlation matrix */ 939 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 940 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 941 for (j=0;j<temp_constraints;j++) { 942 for (k=0;k<j+1;k++) { 943 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)); 944 } 945 } 946 /* compute eigenvalues and eigenvectors of correlation matrix */ 947 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 948 ierr = PetscBLASIntCast(temp_constraints,&Blas_LDA);CHKERRQ(ierr); 949 #if !defined(PETSC_USE_COMPLEX) 950 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,work,&lwork,&lierr)); 951 #else 952 PetscStackCallBLAS("LAPACKsyev",LAPACKsyev_("V","U",&Blas_N,correlation_mat,&Blas_LDA,singular_vals,work,&lwork,rwork,&lierr)); 953 #endif 954 ierr = PetscFPTrapPop();CHKERRQ(ierr); 955 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in SYEV Lapack routine %d",(int)lierr); 956 /* retain eigenvalues greater than tol: note that LAPACKsyev gives eigs in ascending order */ 957 j=0; 958 while (j < temp_constraints && singular_vals[j] < tol) j++; 959 total_counts=total_counts-j; 960 /* scale and copy POD basis into used quadrature memory */ 961 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 962 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 963 ierr = PetscBLASIntCast(temp_constraints,&Blas_K);CHKERRQ(ierr); 964 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 965 ierr = PetscBLASIntCast(temp_constraints,&Blas_LDB);CHKERRQ(ierr); 966 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDC);CHKERRQ(ierr); 967 if (j<temp_constraints) { 968 PetscInt ii; 969 for (k=j;k<temp_constraints;k++) singular_vals[k]=1.0/PetscSqrtReal(singular_vals[k]); 970 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 971 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)); 972 ierr = PetscFPTrapPop();CHKERRQ(ierr); 973 for (k=0;k<temp_constraints-j;k++) { 974 for (ii=0;ii<size_of_constraint;ii++) { 975 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]; 976 } 977 } 978 } 979 #else /* on missing GESVD */ 980 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 981 ierr = PetscBLASIntCast(temp_constraints,&Blas_N);CHKERRQ(ierr); 982 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 983 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 984 #if !defined(PETSC_USE_COMPLEX) 985 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)); 986 #else 987 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)); 988 #endif 989 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GESVD Lapack routine %d",(int)lierr); 990 ierr = PetscFPTrapPop();CHKERRQ(ierr); 991 /* retain eigenvalues greater than tol: note that LAPACKgesvd gives eigs in descending order */ 992 k = temp_constraints; 993 if (k > size_of_constraint) k = size_of_constraint; 994 j = 0; 995 while (j < k && singular_vals[k-j-1] < tol) j++; 996 total_counts = total_counts-temp_constraints+k-j; 997 #endif /* on missing GESVD */ 998 } 999 } 1000 /* free index sets of faces, edges and vertices */ 1001 for (i=0;i<n_ISForFaces;i++) { 1002 ierr = ISDestroy(&ISForFaces[i]);CHKERRQ(ierr); 1003 } 1004 ierr = PetscFree(ISForFaces);CHKERRQ(ierr); 1005 for (i=0;i<n_ISForEdges;i++) { 1006 ierr = ISDestroy(&ISForEdges[i]);CHKERRQ(ierr); 1007 } 1008 ierr = PetscFree(ISForEdges);CHKERRQ(ierr); 1009 ierr = ISDestroy(&ISForVertices);CHKERRQ(ierr); 1010 1011 /* free workspace */ 1012 if (!pcbddc->use_nnsp_true && !skip_lapack) { 1013 ierr = PetscFree(work);CHKERRQ(ierr); 1014 #if defined(PETSC_USE_COMPLEX) 1015 ierr = PetscFree(rwork);CHKERRQ(ierr); 1016 #endif 1017 ierr = PetscFree(singular_vals);CHKERRQ(ierr); 1018 #if defined(PETSC_MISSING_LAPACK_GESVD) 1019 ierr = PetscFree(correlation_mat);CHKERRQ(ierr); 1020 ierr = PetscFree(temp_basis);CHKERRQ(ierr); 1021 #endif 1022 } 1023 for (k=0;k<nnsp_size;k++) { 1024 ierr = VecDestroy(&localnearnullsp[k]);CHKERRQ(ierr); 1025 } 1026 ierr = PetscFree(localnearnullsp);CHKERRQ(ierr); 1027 1028 /* set quantities in pcbddc data structure */ 1029 /* n_vertices defines the number of subdomain corners in the primal space */ 1030 /* n_constraints defines the number of averages (they can be point primal dofs if change of basis is requested) */ 1031 pcbddc->local_primal_size = total_counts; 1032 pcbddc->n_vertices = n_vertices; 1033 pcbddc->n_constraints = pcbddc->local_primal_size-pcbddc->n_vertices; 1034 1035 /* Create constraint matrix */ 1036 /* The constraint matrix is used to compute the l2g map of primal dofs */ 1037 /* so we need to set it up properly either with or without change of basis */ 1038 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ConstraintMatrix);CHKERRQ(ierr); 1039 ierr = MatSetType(pcbddc->ConstraintMatrix,impMatType);CHKERRQ(ierr); 1040 ierr = MatSetSizes(pcbddc->ConstraintMatrix,pcbddc->local_primal_size,pcis->n,pcbddc->local_primal_size,pcis->n);CHKERRQ(ierr); 1041 /* array to compute a local numbering of constraints : vertices first then constraints */ 1042 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&aux_primal_numbering);CHKERRQ(ierr); 1043 /* array to select the proper local node (of minimum index with respect to global ordering) when changing the basis */ 1044 /* 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 */ 1045 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&aux_primal_minloc);CHKERRQ(ierr); 1046 /* auxiliary stuff for basis change */ 1047 ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&global_indices);CHKERRQ(ierr); 1048 ierr = PetscMalloc(pcis->n_B*sizeof(PetscBool),&touched);CHKERRQ(ierr); 1049 ierr = PetscMemzero(touched,pcis->n_B*sizeof(PetscBool));CHKERRQ(ierr); 1050 1051 /* find primal_dofs: subdomain corners plus dofs selected as primal after change of basis */ 1052 total_primal_vertices=0; 1053 for (i=0;i<pcbddc->local_primal_size;i++) { 1054 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 1055 if (size_of_constraint == 1) { 1056 touched[temp_indices_to_constraint_B[temp_indices[i]]]=PETSC_TRUE; 1057 aux_primal_numbering[total_primal_vertices]=temp_indices_to_constraint[temp_indices[i]]; 1058 aux_primal_minloc[total_primal_vertices]=0; 1059 total_primal_vertices++; 1060 } else if (change_basis[i]) { /* Same procedure used in PCBDDCGetPrimalConstraintsLocalIdx */ 1061 PetscInt min_loc,min_index; 1062 ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,&temp_indices_to_constraint[temp_indices[i]],global_indices);CHKERRQ(ierr); 1063 /* find first untouched local node */ 1064 k = 0; 1065 while (touched[temp_indices_to_constraint_B[temp_indices[i]+k]]) k++; 1066 min_index = global_indices[k]; 1067 min_loc = k; 1068 /* search the minimum among global nodes already untouched on the cc */ 1069 for (k=1;k<size_of_constraint;k++) { 1070 /* there can be more than one constraint on a single connected component */ 1071 if (min_index > global_indices[k] && !touched[temp_indices_to_constraint_B[temp_indices[i]+k]]) { 1072 min_index = global_indices[k]; 1073 min_loc = k; 1074 } 1075 } 1076 touched[temp_indices_to_constraint_B[temp_indices[i]+min_loc]] = PETSC_TRUE; 1077 aux_primal_numbering[total_primal_vertices]=temp_indices_to_constraint[temp_indices[i]+min_loc]; 1078 aux_primal_minloc[total_primal_vertices]=min_loc; 1079 total_primal_vertices++; 1080 } 1081 } 1082 /* free workspace */ 1083 ierr = PetscFree(global_indices);CHKERRQ(ierr); 1084 ierr = PetscFree(touched);CHKERRQ(ierr); 1085 /* permute indices in order to have a sorted set of vertices */ 1086 ierr = PetscSortInt(total_primal_vertices,aux_primal_numbering); 1087 1088 /* nonzero structure of constraint matrix */ 1089 ierr = PetscMalloc(pcbddc->local_primal_size*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1090 for (i=0;i<total_primal_vertices;i++) nnz[i]=1; 1091 j=total_primal_vertices; 1092 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1093 if (!change_basis[i]) { 1094 nnz[j]=temp_indices[i+1]-temp_indices[i]; 1095 j++; 1096 } 1097 } 1098 ierr = MatSeqAIJSetPreallocation(pcbddc->ConstraintMatrix,0,nnz);CHKERRQ(ierr); 1099 ierr = PetscFree(nnz);CHKERRQ(ierr); 1100 /* set values in constraint matrix */ 1101 for (i=0;i<total_primal_vertices;i++) { 1102 ierr = MatSetValue(pcbddc->ConstraintMatrix,i,aux_primal_numbering[i],1.0,INSERT_VALUES);CHKERRQ(ierr); 1103 } 1104 total_counts = total_primal_vertices; 1105 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1106 if (!change_basis[i]) { 1107 size_of_constraint=temp_indices[i+1]-temp_indices[i]; 1108 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); 1109 total_counts++; 1110 } 1111 } 1112 /* assembling */ 1113 ierr = MatAssemblyBegin(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1114 ierr = MatAssemblyEnd(pcbddc->ConstraintMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1115 /* 1116 ierr = MatView(pcbddc->ConstraintMatrix,(PetscViewer)0);CHKERRQ(ierr); 1117 */ 1118 /* Create matrix for change of basis. We don't need it in case pcbddc->use_change_of_basis is FALSE */ 1119 if (pcbddc->use_change_of_basis) { 1120 PetscBool qr_needed = PETSC_FALSE; 1121 /* change of basis acts on local interfaces -> dimension is n_B x n_B */ 1122 ierr = MatCreate(PETSC_COMM_SELF,&pcbddc->ChangeOfBasisMatrix);CHKERRQ(ierr); 1123 ierr = MatSetType(pcbddc->ChangeOfBasisMatrix,impMatType);CHKERRQ(ierr); 1124 ierr = MatSetSizes(pcbddc->ChangeOfBasisMatrix,pcis->n_B,pcis->n_B,pcis->n_B,pcis->n_B);CHKERRQ(ierr); 1125 /* work arrays */ 1126 ierr = PetscMalloc(pcis->n_B*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1127 for (i=0;i<pcis->n_B;i++) nnz[i]=1; 1128 /* nonzeros per row */ 1129 for (i=pcbddc->n_vertices;i<pcbddc->local_primal_size;i++) { 1130 if (change_basis[i]) { 1131 qr_needed = PETSC_TRUE; 1132 size_of_constraint = temp_indices[i+1]-temp_indices[i]; 1133 for (j=0;j<size_of_constraint;j++) nnz[temp_indices_to_constraint_B[temp_indices[i]+j]] = size_of_constraint; 1134 } 1135 } 1136 ierr = MatSeqAIJSetPreallocation(pcbddc->ChangeOfBasisMatrix,0,nnz);CHKERRQ(ierr); 1137 ierr = PetscFree(nnz);CHKERRQ(ierr); 1138 /* Set initial identity in the matrix */ 1139 for (i=0;i<pcis->n_B;i++) { 1140 ierr = MatSetValue(pcbddc->ChangeOfBasisMatrix,i,i,1.0,INSERT_VALUES);CHKERRQ(ierr); 1141 } 1142 1143 /* Now we loop on the constraints which need a change of basis */ 1144 /* Change of basis matrix is evaluated as the FIRST APPROACH in */ 1145 /* Klawonn and Widlund, Dual-primal FETI-DP methods for linear elasticity, (see Sect 6.2.1) */ 1146 /* Change of basis matrix T computed via QR decomposition of constraints */ 1147 if (qr_needed) { 1148 /* dual and primal dofs on a single cc */ 1149 PetscInt dual_dofs,primal_dofs; 1150 /* iterator on aux_primal_minloc (ordered as read from nearnullspace: vertices, edges and then constraints) */ 1151 PetscInt primal_counter; 1152 /* working stuff for GEQRF */ 1153 PetscScalar *qr_basis,*qr_tau,*qr_work,lqr_work_t; 1154 PetscBLASInt lqr_work; 1155 /* working stuff for UNGQR */ 1156 PetscScalar *gqr_work,lgqr_work_t; 1157 PetscBLASInt lgqr_work; 1158 /* working stuff for TRTRS */ 1159 PetscScalar *trs_rhs; 1160 PetscBLASInt Blas_NRHS; 1161 /* pointers for values insertion into change of basis matrix */ 1162 PetscInt *start_rows,*start_cols; 1163 PetscScalar *start_vals; 1164 /* working stuff for values insertion */ 1165 PetscBool *is_primal; 1166 1167 /* space to store Q */ 1168 ierr = PetscMalloc((max_size_of_constraint)*(max_size_of_constraint)*sizeof(PetscScalar),&qr_basis);CHKERRQ(ierr); 1169 /* first we issue queries for optimal work */ 1170 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_M);CHKERRQ(ierr); 1171 ierr = PetscBLASIntCast(max_constraints,&Blas_N);CHKERRQ(ierr); 1172 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1173 lqr_work = -1; 1174 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,&lqr_work_t,&lqr_work,&lierr)); 1175 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to GEQRF Lapack routine %d",(int)lierr); 1176 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lqr_work_t),&lqr_work);CHKERRQ(ierr); 1177 ierr = PetscMalloc((PetscInt)PetscRealPart(lqr_work_t)*sizeof(*qr_work),&qr_work);CHKERRQ(ierr); 1178 lgqr_work = -1; 1179 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_M);CHKERRQ(ierr); 1180 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_N);CHKERRQ(ierr); 1181 ierr = PetscBLASIntCast(max_constraints,&Blas_K);CHKERRQ(ierr); 1182 ierr = PetscBLASIntCast(max_size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1183 if (Blas_K>Blas_M) Blas_K=Blas_M; /* adjust just for computing optimal work */ 1184 PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Blas_M,&Blas_N,&Blas_K,qr_basis,&Blas_LDA,qr_tau,&lgqr_work_t,&lgqr_work,&lierr)); 1185 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in query to UNGQR Lapack routine %d",(int)lierr); 1186 ierr = PetscBLASIntCast((PetscInt)PetscRealPart(lgqr_work_t),&lgqr_work);CHKERRQ(ierr); 1187 ierr = PetscMalloc((PetscInt)PetscRealPart(lgqr_work_t)*sizeof(*gqr_work),&gqr_work);CHKERRQ(ierr); 1188 /* array to store scaling factors for reflectors */ 1189 ierr = PetscMalloc(max_constraints*sizeof(*qr_tau),&qr_tau);CHKERRQ(ierr); 1190 /* array to store rhs and solution of triangular solver */ 1191 ierr = PetscMalloc(max_constraints*max_constraints*sizeof(*trs_rhs),&trs_rhs);CHKERRQ(ierr); 1192 /* array to store whether a node is primal or not */ 1193 ierr = PetscMalloc(pcis->n_B*sizeof(*is_primal),&is_primal);CHKERRQ(ierr); 1194 ierr = PetscMemzero(is_primal,pcis->n_B*sizeof(*is_primal));CHKERRQ(ierr); 1195 for (i=0;i<total_primal_vertices;i++) is_primal[local_to_B[aux_primal_numbering[i]]] = PETSC_TRUE; 1196 1197 /* allocating workspace for check */ 1198 if (pcbddc->dbg_flag) { 1199 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"--------------------------------------------------------------\n");CHKERRQ(ierr); 1200 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"Checking change of basis computation for subdomain %04d\n",PetscGlobalRank);CHKERRQ(ierr); 1201 ierr = PetscMalloc(max_size_of_constraint*(max_constraints+max_size_of_constraint)*sizeof(*work),&work);CHKERRQ(ierr); 1202 } 1203 1204 /* loop on constraints and see whether or not they need a change of basis */ 1205 /* -> using implicit ordering contained in temp_indices data */ 1206 total_counts = pcbddc->n_vertices; 1207 primal_counter = total_counts; 1208 while (total_counts<pcbddc->local_primal_size) { 1209 primal_dofs = 1; 1210 if (change_basis[total_counts]) { 1211 /* get all constraints with same support: if more then one constraint is present on the cc then surely indices are stored contiguosly */ 1212 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]]) { 1213 primal_dofs++; 1214 } 1215 /* get constraint info */ 1216 size_of_constraint = temp_indices[total_counts+1]-temp_indices[total_counts]; 1217 dual_dofs = size_of_constraint-primal_dofs; 1218 1219 /* copy quadrature constraints for change of basis check */ 1220 if (pcbddc->dbg_flag) { 1221 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); 1222 ierr = PetscMemcpy(work,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1223 } 1224 1225 /* copy temporary constraints into larger work vector (in order to store all columns of Q) */ 1226 ierr = PetscMemcpy(qr_basis,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1227 1228 /* compute QR decomposition of constraints */ 1229 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1230 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1231 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1232 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1233 PetscStackCallBLAS("LAPACKgeqrf",LAPACKgeqrf_(&Blas_M,&Blas_N,qr_basis,&Blas_LDA,qr_tau,qr_work,&lqr_work,&lierr)); 1234 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in GEQRF Lapack routine %d",(int)lierr); 1235 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1236 1237 /* explictly compute R^-T */ 1238 ierr = PetscMemzero(trs_rhs,primal_dofs*primal_dofs*sizeof(*trs_rhs));CHKERRQ(ierr); 1239 for (j=0;j<primal_dofs;j++) trs_rhs[j*(primal_dofs+1)] = 1.0; 1240 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1241 ierr = PetscBLASIntCast(primal_dofs,&Blas_NRHS);CHKERRQ(ierr); 1242 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1243 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDB);CHKERRQ(ierr); 1244 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1245 PetscStackCallBLAS("LAPACKtrtrs",LAPACKtrtrs_("U","T","N",&Blas_N,&Blas_NRHS,qr_basis,&Blas_LDA,trs_rhs,&Blas_LDB,&lierr)); 1246 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in TRTRS Lapack routine %d",(int)lierr); 1247 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1248 1249 /* explcitly compute all columns of Q (Q = [Q1 | Q2] ) overwriting QR factorization in qr_basis */ 1250 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1251 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1252 ierr = PetscBLASIntCast(primal_dofs,&Blas_K);CHKERRQ(ierr); 1253 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1254 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1255 PetscStackCallBLAS("LAPACKungqr",LAPACKungqr_(&Blas_M,&Blas_N,&Blas_K,qr_basis,&Blas_LDA,qr_tau,gqr_work,&lgqr_work,&lierr)); 1256 if (lierr) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in UNGQR Lapack routine %d",(int)lierr); 1257 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1258 1259 /* first primal_dofs columns of Q need to be re-scaled in order to be unitary w.r.t constraints 1260 i.e. C_{pxn}*Q_{nxn} should be equal to [I_pxp | 0_pxd] (see check below) 1261 where n=size_of_constraint, p=primal_dofs, d=dual_dofs (n=p+d), I and 0 identity and null matrix resp. */ 1262 ierr = PetscBLASIntCast(size_of_constraint,&Blas_M);CHKERRQ(ierr); 1263 ierr = PetscBLASIntCast(primal_dofs,&Blas_N);CHKERRQ(ierr); 1264 ierr = PetscBLASIntCast(primal_dofs,&Blas_K);CHKERRQ(ierr); 1265 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1266 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDB);CHKERRQ(ierr); 1267 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDC);CHKERRQ(ierr); 1268 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1269 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)); 1270 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1271 ierr = PetscMemcpy(qr_basis,&temp_quadrature_constraint[temp_indices[total_counts]],size_of_constraint*primal_dofs*sizeof(PetscScalar));CHKERRQ(ierr); 1272 1273 /* insert values in change of basis matrix respecting global ordering of new primal dofs */ 1274 start_rows = &temp_indices_to_constraint_B[temp_indices[total_counts]]; 1275 /* insert cols for primal dofs */ 1276 for (j=0;j<primal_dofs;j++) { 1277 start_vals = &qr_basis[j*size_of_constraint]; 1278 start_cols = &temp_indices_to_constraint_B[temp_indices[total_counts]+aux_primal_minloc[primal_counter+j]]; 1279 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,start_rows,1,start_cols,start_vals,INSERT_VALUES);CHKERRQ(ierr); 1280 } 1281 /* insert cols for dual dofs */ 1282 for (j=0,k=0;j<dual_dofs;k++) { 1283 if (!is_primal[temp_indices_to_constraint_B[temp_indices[total_counts]+k]]) { 1284 start_vals = &qr_basis[(primal_dofs+j)*size_of_constraint]; 1285 start_cols = &temp_indices_to_constraint_B[temp_indices[total_counts]+k]; 1286 ierr = MatSetValues(pcbddc->ChangeOfBasisMatrix,size_of_constraint,start_rows,1,start_cols,start_vals,INSERT_VALUES);CHKERRQ(ierr); 1287 j++; 1288 } 1289 } 1290 1291 /* check change of basis */ 1292 if (pcbddc->dbg_flag) { 1293 PetscInt ii,jj; 1294 PetscBool valid_qr=PETSC_TRUE; 1295 ierr = PetscBLASIntCast(primal_dofs,&Blas_M);CHKERRQ(ierr); 1296 ierr = PetscBLASIntCast(size_of_constraint,&Blas_N);CHKERRQ(ierr); 1297 ierr = PetscBLASIntCast(size_of_constraint,&Blas_K);CHKERRQ(ierr); 1298 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDA);CHKERRQ(ierr); 1299 ierr = PetscBLASIntCast(size_of_constraint,&Blas_LDB);CHKERRQ(ierr); 1300 ierr = PetscBLASIntCast(primal_dofs,&Blas_LDC);CHKERRQ(ierr); 1301 ierr = PetscFPTrapPush(PETSC_FP_TRAP_OFF);CHKERRQ(ierr); 1302 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)); 1303 ierr = PetscFPTrapPop();CHKERRQ(ierr); 1304 for (jj=0;jj<size_of_constraint;jj++) { 1305 for (ii=0;ii<primal_dofs;ii++) { 1306 if (ii != jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]) > 1.e-12) valid_qr = PETSC_FALSE; 1307 if (ii == jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]-1.0) > 1.e-12) valid_qr = PETSC_FALSE; 1308 } 1309 } 1310 if (!valid_qr) { 1311 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\t-> wrong change of basis!\n",PetscGlobalRank);CHKERRQ(ierr); 1312 for (jj=0;jj<size_of_constraint;jj++) { 1313 for (ii=0;ii<primal_dofs;ii++) { 1314 if (ii != jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]) > 1.e-12) { 1315 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])); 1316 } 1317 if (ii == jj && PetscAbsScalar(work[size_of_constraint*primal_dofs+jj*primal_dofs+ii]-1.0) > 1.e-12) { 1318 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])); 1319 } 1320 } 1321 } 1322 } else { 1323 ierr = PetscViewerASCIISynchronizedPrintf(pcbddc->dbg_viewer,"\t-> right change of basis!\n",PetscGlobalRank);CHKERRQ(ierr); 1324 } 1325 } 1326 /* increment primal counter */ 1327 primal_counter += primal_dofs; 1328 } else { 1329 if (pcbddc->dbg_flag) { 1330 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); 1331 } 1332 } 1333 /* increment constraint counter total_counts */ 1334 total_counts += primal_dofs; 1335 } 1336 if (pcbddc->dbg_flag) { 1337 ierr = PetscViewerFlush(pcbddc->dbg_viewer);CHKERRQ(ierr); 1338 ierr = PetscFree(work);CHKERRQ(ierr); 1339 } 1340 /* free workspace */ 1341 ierr = PetscFree(trs_rhs);CHKERRQ(ierr); 1342 ierr = PetscFree(qr_tau);CHKERRQ(ierr); 1343 ierr = PetscFree(qr_work);CHKERRQ(ierr); 1344 ierr = PetscFree(gqr_work);CHKERRQ(ierr); 1345 ierr = PetscFree(is_primal);CHKERRQ(ierr); 1346 ierr = PetscFree(qr_basis);CHKERRQ(ierr); 1347 } 1348 /* assembling */ 1349 ierr = MatAssemblyBegin(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1350 ierr = MatAssemblyEnd(pcbddc->ChangeOfBasisMatrix,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1351 /* 1352 ierr = MatView(pcbddc->ChangeOfBasisMatrix,(PetscViewer)0);CHKERRQ(ierr); 1353 */ 1354 } 1355 /* free workspace */ 1356 ierr = PetscFree(aux_primal_numbering);CHKERRQ(ierr); 1357 ierr = PetscFree(aux_primal_minloc);CHKERRQ(ierr); 1358 ierr = PetscFree(temp_indices);CHKERRQ(ierr); 1359 ierr = PetscFree(change_basis);CHKERRQ(ierr); 1360 ierr = PetscFree(temp_indices_to_constraint);CHKERRQ(ierr); 1361 ierr = PetscFree(temp_indices_to_constraint_B);CHKERRQ(ierr); 1362 ierr = PetscFree(local_to_B);CHKERRQ(ierr); 1363 ierr = PetscFree(temp_quadrature_constraint);CHKERRQ(ierr); 1364 PetscFunctionReturn(0); 1365 } 1366 1367 #undef __FUNCT__ 1368 #define __FUNCT__ "PCBDDCAnalyzeInterface" 1369 PetscErrorCode PCBDDCAnalyzeInterface(PC pc) 1370 { 1371 PC_BDDC *pcbddc = (PC_BDDC*)pc->data; 1372 PC_IS *pcis = (PC_IS*)pc->data; 1373 Mat_IS *matis = (Mat_IS*)pc->pmat->data; 1374 PetscInt bs,ierr,i,vertex_size; 1375 PetscViewer viewer=pcbddc->dbg_viewer; 1376 1377 PetscFunctionBegin; 1378 /* Init local Graph struct */ 1379 ierr = PCBDDCGraphInit(pcbddc->mat_graph,matis->mapping);CHKERRQ(ierr); 1380 1381 /* Check validity of the csr graph passed in by the user */ 1382 if (pcbddc->mat_graph->nvtxs_csr != pcbddc->mat_graph->nvtxs) { 1383 ierr = PCBDDCGraphResetCSR(pcbddc->mat_graph);CHKERRQ(ierr); 1384 } 1385 /* Set default CSR adjacency of local dofs if not provided by the user with PCBDDCSetLocalAdjacencyGraph */ 1386 if (!pcbddc->mat_graph->xadj || !pcbddc->mat_graph->adjncy) { 1387 Mat mat_adj; 1388 const PetscInt *xadj,*adjncy; 1389 PetscBool flg_row=PETSC_TRUE; 1390 1391 ierr = MatConvert(matis->A,MATMPIADJ,MAT_INITIAL_MATRIX,&mat_adj);CHKERRQ(ierr); 1392 ierr = MatGetRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1393 if (!flg_row) { 1394 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatGetRowIJ called in %s\n",__FUNCT__); 1395 } 1396 ierr = PCBDDCSetLocalAdjacencyGraph(pc,i,xadj,adjncy,PETSC_COPY_VALUES);CHKERRQ(ierr); 1397 ierr = MatRestoreRowIJ(mat_adj,0,PETSC_TRUE,PETSC_FALSE,&i,&xadj,&adjncy,&flg_row);CHKERRQ(ierr); 1398 if (!flg_row) { 1399 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Error in MatRestoreRowIJ called in %s\n",__FUNCT__); 1400 } 1401 ierr = MatDestroy(&mat_adj);CHKERRQ(ierr); 1402 } 1403 1404 /* Set default dofs' splitting if no information has been provided by the user with PCBDDCSetDofsSplitting */ 1405 vertex_size = 1; 1406 if (!pcbddc->n_ISForDofs) { 1407 IS *custom_ISForDofs; 1408 1409 ierr = MatGetBlockSize(matis->A,&bs);CHKERRQ(ierr); 1410 ierr = PetscMalloc(bs*sizeof(IS),&custom_ISForDofs);CHKERRQ(ierr); 1411 for (i=0;i<bs;i++) { 1412 ierr = ISCreateStride(PETSC_COMM_SELF,pcis->n/bs,i,bs,&custom_ISForDofs[i]);CHKERRQ(ierr); 1413 } 1414 ierr = PCBDDCSetDofsSplitting(pc,bs,custom_ISForDofs);CHKERRQ(ierr); 1415 /* remove my references to IS objects */ 1416 for (i=0;i<bs;i++) { 1417 ierr = ISDestroy(&custom_ISForDofs[i]);CHKERRQ(ierr); 1418 } 1419 ierr = PetscFree(custom_ISForDofs);CHKERRQ(ierr); 1420 } else { /* mat block size as vertex size (used for elasticity) */ 1421 ierr = MatGetBlockSize(matis->A,&vertex_size);CHKERRQ(ierr); 1422 } 1423 1424 /* Setup of Graph */ 1425 ierr = PCBDDCGraphSetUp(pcbddc->mat_graph,vertex_size,pcbddc->NeumannBoundaries,pcbddc->DirichletBoundaries,pcbddc->n_ISForDofs,pcbddc->ISForDofs,pcbddc->user_primal_vertices); 1426 1427 /* Graph's connected components analysis */ 1428 ierr = PCBDDCGraphComputeConnectedComponents(pcbddc->mat_graph);CHKERRQ(ierr); 1429 1430 /* print some info to stdout */ 1431 if (pcbddc->dbg_flag) { 1432 ierr = PCBDDCGraphASCIIView(pcbddc->mat_graph,pcbddc->dbg_flag,viewer); 1433 } 1434 PetscFunctionReturn(0); 1435 } 1436 1437 #undef __FUNCT__ 1438 #define __FUNCT__ "PCBDDCGetPrimalVerticesLocalIdx" 1439 PetscErrorCode PCBDDCGetPrimalVerticesLocalIdx(PC pc, PetscInt *n_vertices, PetscInt *vertices_idx[]) 1440 { 1441 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1442 PetscInt *vertices,*row_cmat_indices,n,i,size_of_constraint,local_primal_size; 1443 PetscErrorCode ierr; 1444 1445 PetscFunctionBegin; 1446 n = 0; 1447 vertices = 0; 1448 if (pcbddc->ConstraintMatrix) { 1449 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&i);CHKERRQ(ierr); 1450 for (i=0;i<local_primal_size;i++) { 1451 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1452 if (size_of_constraint == 1) n++; 1453 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1454 } 1455 if (vertices_idx) { 1456 ierr = PetscMalloc(n*sizeof(PetscInt),&vertices);CHKERRQ(ierr); 1457 n = 0; 1458 for (i=0;i<local_primal_size;i++) { 1459 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1460 if (size_of_constraint == 1) { 1461 vertices[n++]=row_cmat_indices[0]; 1462 } 1463 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1464 } 1465 } 1466 } 1467 *n_vertices = n; 1468 if (vertices_idx) *vertices_idx = vertices; 1469 PetscFunctionReturn(0); 1470 } 1471 1472 #undef __FUNCT__ 1473 #define __FUNCT__ "PCBDDCGetPrimalConstraintsLocalIdx" 1474 PetscErrorCode PCBDDCGetPrimalConstraintsLocalIdx(PC pc, PetscInt *n_constraints, PetscInt *constraints_idx[]) 1475 { 1476 PC_BDDC *pcbddc = (PC_BDDC*)(pc->data); 1477 PetscInt *constraints_index,*row_cmat_indices,*row_cmat_global_indices; 1478 PetscInt n,i,j,size_of_constraint,local_primal_size,local_size,max_size_of_constraint,min_index,min_loc; 1479 PetscBool *touched; 1480 PetscErrorCode ierr; 1481 1482 PetscFunctionBegin; 1483 n = 0; 1484 constraints_index = 0; 1485 if (pcbddc->ConstraintMatrix) { 1486 ierr = MatGetSize(pcbddc->ConstraintMatrix,&local_primal_size,&local_size);CHKERRQ(ierr); 1487 max_size_of_constraint = 0; 1488 for (i=0;i<local_primal_size;i++) { 1489 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1490 if (size_of_constraint > 1) { 1491 n++; 1492 } 1493 max_size_of_constraint = PetscMax(size_of_constraint,max_size_of_constraint); 1494 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,NULL,NULL);CHKERRQ(ierr); 1495 } 1496 if (constraints_idx) { 1497 ierr = PetscMalloc(n*sizeof(PetscInt),&constraints_index);CHKERRQ(ierr); 1498 ierr = PetscMalloc(max_size_of_constraint*sizeof(PetscInt),&row_cmat_global_indices);CHKERRQ(ierr); 1499 ierr = PetscMalloc(local_size*sizeof(PetscBool),&touched);CHKERRQ(ierr); 1500 ierr = PetscMemzero(touched,local_size*sizeof(PetscBool));CHKERRQ(ierr); 1501 n = 0; 1502 for (i=0;i<local_primal_size;i++) { 1503 ierr = MatGetRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1504 if (size_of_constraint > 1) { 1505 ierr = ISLocalToGlobalMappingApply(pcbddc->mat_graph->l2gmap,size_of_constraint,row_cmat_indices,row_cmat_global_indices);CHKERRQ(ierr); 1506 /* find first untouched local node */ 1507 j = 0; 1508 while(touched[row_cmat_indices[j]]) j++; 1509 min_index = row_cmat_global_indices[j]; 1510 min_loc = j; 1511 /* search the minimum among nodes not yet touched on the connected component 1512 since there can be more than one constraint on a single cc */ 1513 for (j=1;j<size_of_constraint;j++) { 1514 if (min_index > row_cmat_global_indices[j] && !touched[row_cmat_indices[j]]) { 1515 min_index = row_cmat_global_indices[j]; 1516 min_loc = j; 1517 } 1518 } 1519 touched[row_cmat_indices[min_loc]] = PETSC_TRUE; 1520 constraints_index[n++] = row_cmat_indices[min_loc]; 1521 } 1522 ierr = MatRestoreRow(pcbddc->ConstraintMatrix,i,&size_of_constraint,(const PetscInt**)&row_cmat_indices,NULL);CHKERRQ(ierr); 1523 } 1524 ierr = PetscFree(touched);CHKERRQ(ierr); 1525 ierr = PetscFree(row_cmat_global_indices);CHKERRQ(ierr); 1526 } 1527 } 1528 *n_constraints = n; 1529 if (constraints_idx) *constraints_idx = constraints_index; 1530 PetscFunctionReturn(0); 1531 } 1532 1533 /* the next two functions has been adapted from pcis.c */ 1534 #undef __FUNCT__ 1535 #define __FUNCT__ "PCBDDCApplySchur" 1536 PetscErrorCode PCBDDCApplySchur(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 1537 { 1538 PetscErrorCode ierr; 1539 PC_IS *pcis = (PC_IS*)(pc->data); 1540 1541 PetscFunctionBegin; 1542 if (!vec2_B) { vec2_B = v; } 1543 ierr = MatMult(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 1544 ierr = MatMult(pcis->A_IB,v,vec1_D);CHKERRQ(ierr); 1545 ierr = KSPSolve(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 1546 ierr = MatMult(pcis->A_BI,vec2_D,vec2_B);CHKERRQ(ierr); 1547 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 1548 PetscFunctionReturn(0); 1549 } 1550 1551 #undef __FUNCT__ 1552 #define __FUNCT__ "PCBDDCApplySchurTranspose" 1553 PetscErrorCode PCBDDCApplySchurTranspose(PC pc, Vec v, Vec vec1_B, Vec vec2_B, Vec vec1_D, Vec vec2_D) 1554 { 1555 PetscErrorCode ierr; 1556 PC_IS *pcis = (PC_IS*)(pc->data); 1557 1558 PetscFunctionBegin; 1559 if (!vec2_B) { vec2_B = v; } 1560 ierr = MatMultTranspose(pcis->A_BB,v,vec1_B);CHKERRQ(ierr); 1561 ierr = MatMultTranspose(pcis->A_BI,v,vec1_D);CHKERRQ(ierr); 1562 ierr = KSPSolveTranspose(pcis->ksp_D,vec1_D,vec2_D);CHKERRQ(ierr); 1563 ierr = MatMultTranspose(pcis->A_IB,vec2_D,vec2_B);CHKERRQ(ierr); 1564 ierr = VecAXPY(vec1_B,-1.0,vec2_B);CHKERRQ(ierr); 1565 PetscFunctionReturn(0); 1566 } 1567 1568 #undef __FUNCT__ 1569 #define __FUNCT__ "PCBDDCSubsetNumbering" 1570 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[]) 1571 { 1572 Vec local_vec,global_vec; 1573 IS seqis,paris; 1574 VecScatter scatter_ctx; 1575 PetscScalar *array; 1576 PetscInt *temp_global_dofs; 1577 PetscScalar globalsum; 1578 PetscInt i,j,s; 1579 PetscInt nlocals,first_index,old_index,max_local; 1580 PetscMPIInt rank_prec_comm,size_prec_comm,max_global; 1581 PetscMPIInt *dof_sizes,*dof_displs; 1582 PetscBool first_found; 1583 PetscErrorCode ierr; 1584 1585 PetscFunctionBegin; 1586 /* mpi buffers */ 1587 MPI_Comm_size(comm,&size_prec_comm); 1588 MPI_Comm_rank(comm,&rank_prec_comm); 1589 j = ( !rank_prec_comm ? size_prec_comm : 0); 1590 ierr = PetscMalloc(j*sizeof(*dof_sizes),&dof_sizes);CHKERRQ(ierr); 1591 ierr = PetscMalloc(j*sizeof(*dof_displs),&dof_displs);CHKERRQ(ierr); 1592 /* get maximum size of subset */ 1593 ierr = PetscMalloc(n_local_dofs*sizeof(PetscInt),&temp_global_dofs);CHKERRQ(ierr); 1594 ierr = ISLocalToGlobalMappingApply(l2gmap,n_local_dofs,local_dofs,temp_global_dofs);CHKERRQ(ierr); 1595 max_local = 0; 1596 if (n_local_dofs) { 1597 max_local = temp_global_dofs[0]; 1598 for (i=1;i<n_local_dofs;i++) { 1599 if (max_local < temp_global_dofs[i] ) { 1600 max_local = temp_global_dofs[i]; 1601 } 1602 } 1603 } 1604 ierr = MPI_Allreduce(&max_local,&max_global,1,MPIU_INT,MPI_MAX,comm); 1605 max_global++; 1606 max_local = 0; 1607 if (n_local_dofs) { 1608 max_local = local_dofs[0]; 1609 for (i=1;i<n_local_dofs;i++) { 1610 if (max_local < local_dofs[i] ) { 1611 max_local = local_dofs[i]; 1612 } 1613 } 1614 } 1615 max_local++; 1616 /* allocate workspace */ 1617 ierr = VecCreate(PETSC_COMM_SELF,&local_vec);CHKERRQ(ierr); 1618 ierr = VecSetSizes(local_vec,PETSC_DECIDE,max_local);CHKERRQ(ierr); 1619 ierr = VecSetType(local_vec,VECSEQ);CHKERRQ(ierr); 1620 ierr = VecCreate(comm,&global_vec);CHKERRQ(ierr); 1621 ierr = VecSetSizes(global_vec,PETSC_DECIDE,max_global);CHKERRQ(ierr); 1622 ierr = VecSetType(global_vec,VECMPI);CHKERRQ(ierr); 1623 /* create scatter */ 1624 ierr = ISCreateGeneral(PETSC_COMM_SELF,n_local_dofs,local_dofs,PETSC_COPY_VALUES,&seqis);CHKERRQ(ierr); 1625 ierr = ISCreateGeneral(comm,n_local_dofs,temp_global_dofs,PETSC_COPY_VALUES,&paris);CHKERRQ(ierr); 1626 ierr = VecScatterCreate(local_vec,seqis,global_vec,paris,&scatter_ctx);CHKERRQ(ierr); 1627 ierr = ISDestroy(&seqis);CHKERRQ(ierr); 1628 ierr = ISDestroy(&paris);CHKERRQ(ierr); 1629 /* init array */ 1630 ierr = VecSet(global_vec,0.0);CHKERRQ(ierr); 1631 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 1632 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 1633 if (local_dofs_mult) { 1634 for (i=0;i<n_local_dofs;i++) { 1635 array[local_dofs[i]]=(PetscScalar)local_dofs_mult[i]; 1636 } 1637 } else { 1638 for (i=0;i<n_local_dofs;i++) { 1639 array[local_dofs[i]]=1.0; 1640 } 1641 } 1642 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 1643 /* scatter into global vec and get total number of global dofs */ 1644 ierr = VecScatterBegin(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1645 ierr = VecScatterEnd(scatter_ctx,local_vec,global_vec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 1646 ierr = VecSum(global_vec,&globalsum);CHKERRQ(ierr); 1647 *n_global_subset = (PetscInt)PetscRealPart(globalsum); 1648 /* Fill global_vec with cumulative function for global numbering */ 1649 ierr = VecGetArray(global_vec,&array);CHKERRQ(ierr); 1650 ierr = VecGetLocalSize(global_vec,&s);CHKERRQ(ierr); 1651 nlocals = 0; 1652 first_index = -1; 1653 first_found = PETSC_FALSE; 1654 for (i=0;i<s;i++) { 1655 if (!first_found && PetscRealPart(array[i]) > 0.0) { 1656 first_found = PETSC_TRUE; 1657 first_index = i; 1658 } 1659 nlocals += (PetscInt)PetscRealPart(array[i]); 1660 } 1661 ierr = MPI_Gather(&nlocals,1,MPIU_INT,dof_sizes,1,MPIU_INT,0,comm);CHKERRQ(ierr); 1662 if (!rank_prec_comm) { 1663 dof_displs[0]=0; 1664 for (i=1;i<size_prec_comm;i++) { 1665 dof_displs[i] = dof_displs[i-1]+dof_sizes[i-1]; 1666 } 1667 } 1668 ierr = MPI_Scatter(dof_displs,1,MPIU_INT,&nlocals,1,MPIU_INT,0,comm);CHKERRQ(ierr); 1669 if (first_found) { 1670 array[first_index] += (PetscScalar)nlocals; 1671 old_index = first_index; 1672 for (i=first_index+1;i<s;i++) { 1673 if (PetscRealPart(array[i]) > 0.0) { 1674 array[i] += array[old_index]; 1675 old_index = i; 1676 } 1677 } 1678 } 1679 ierr = VecRestoreArray(global_vec,&array);CHKERRQ(ierr); 1680 ierr = VecSet(local_vec,0.0);CHKERRQ(ierr); 1681 ierr = VecScatterBegin(scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1682 ierr = VecScatterEnd (scatter_ctx,global_vec,local_vec,INSERT_VALUES,SCATTER_REVERSE);CHKERRQ(ierr); 1683 /* get global ordering of local dofs */ 1684 ierr = VecGetArray(local_vec,&array);CHKERRQ(ierr); 1685 if (local_dofs_mult) { 1686 for (i=0;i<n_local_dofs;i++) { 1687 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-local_dofs_mult[i]; 1688 } 1689 } else { 1690 for (i=0;i<n_local_dofs;i++) { 1691 temp_global_dofs[i] = (PetscInt)PetscRealPart(array[local_dofs[i]])-1; 1692 } 1693 } 1694 ierr = VecRestoreArray(local_vec,&array);CHKERRQ(ierr); 1695 /* free workspace */ 1696 ierr = VecScatterDestroy(&scatter_ctx);CHKERRQ(ierr); 1697 ierr = VecDestroy(&local_vec);CHKERRQ(ierr); 1698 ierr = VecDestroy(&global_vec);CHKERRQ(ierr); 1699 ierr = PetscFree(dof_sizes);CHKERRQ(ierr); 1700 ierr = PetscFree(dof_displs);CHKERRQ(ierr); 1701 /* return pointer to global ordering of local dofs */ 1702 *global_numbering_subset = temp_global_dofs; 1703 PetscFunctionReturn(0); 1704 } 1705 1706 #undef __FUNCT__ 1707 #define __FUNCT__ "PCBDDCOrthonormalizeVecs" 1708 PetscErrorCode PCBDDCOrthonormalizeVecs(PetscInt n, Vec vecs[]) 1709 { 1710 PetscInt i,j; 1711 PetscScalar *alphas; 1712 PetscErrorCode ierr; 1713 1714 PetscFunctionBegin; 1715 /* this implements stabilized Gram-Schmidt */ 1716 ierr = PetscMalloc(n*sizeof(PetscScalar),&alphas);CHKERRQ(ierr); 1717 for (i=0;i<n;i++) { 1718 ierr = VecNormalize(vecs[i],NULL);CHKERRQ(ierr); 1719 if (i<n) { ierr = VecMDot(vecs[i],n-i-1,&vecs[i+1],&alphas[i+1]);CHKERRQ(ierr); } 1720 for (j=i+1;j<n;j++) { ierr = VecAXPY(vecs[j],PetscConj(-alphas[j]),vecs[i]);CHKERRQ(ierr); } 1721 } 1722 ierr = PetscFree(alphas);CHKERRQ(ierr); 1723 PetscFunctionReturn(0); 1724 } 1725 1726 1727