1 /* 2 GAMG geometric-algebric multigrid PC - Mark Adams 2011 3 */ 4 #include <petsc/private/matimpl.h> 5 #include <../src/ksp/pc/impls/gamg/gamg.h> /*I "petscpc.h" I*/ 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatCollapseRow" 9 /* 10 Produces a set of block column indices of the matrix row, one for each block represented in the original row 11 12 n - the number of block indices in cc[] 13 cc - the block indices (must be large enough to contain the indices) 14 */ 15 PETSC_STATIC_INLINE PetscErrorCode MatCollapseRow(Mat Amat,PetscInt row,PetscInt bs,PetscInt *n,PetscInt *cc) 16 { 17 PetscInt cnt = -1,nidx,j; 18 const PetscInt *idx; 19 PetscErrorCode ierr; 20 21 PetscFunctionBegin; 22 ierr = MatGetRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr); 23 if (nidx) { 24 cnt = 0; 25 cc[cnt] = idx[0]/bs; 26 for (j=1; j<nidx; j++) { 27 if (cc[cnt] < idx[j]/bs) cc[++cnt] = idx[j]/bs; 28 } 29 } 30 ierr = MatRestoreRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr); 31 *n = cnt+1; 32 PetscFunctionReturn(0); 33 } 34 35 #undef __FUNCT__ 36 #define __FUNCT__ "MatCollapseRows" 37 /* 38 Produces a set of block column indices of the matrix block row, one for each block represented in the original set of rows 39 40 ncollapsed - the number of block indices 41 collapsed - the block indices (must be large enough to contain the indices) 42 */ 43 PETSC_STATIC_INLINE PetscErrorCode MatCollapseRows(Mat Amat,PetscInt start,PetscInt bs,PetscInt *w0,PetscInt *w1,PetscInt *w2,PetscInt *ncollapsed,PetscInt **collapsed) 44 { 45 PetscInt i,nprev,*cprev = w0,ncur = 0,*ccur = w1,*merged = w2,*cprevtmp; 46 PetscErrorCode ierr; 47 48 PetscFunctionBegin; 49 ierr = MatCollapseRow(Amat,start,bs,&nprev,cprev);CHKERRQ(ierr); 50 for (i=start+1; i<start+bs; i++) { 51 ierr = MatCollapseRow(Amat,i,bs,&ncur,ccur);CHKERRQ(ierr); 52 ierr = PetscMergeIntArray(nprev,cprev,ncur,ccur,&nprev,&merged);CHKERRQ(ierr); 53 cprevtmp = cprev; cprev = merged; merged = cprevtmp; 54 } 55 *ncollapsed = nprev; 56 if (collapsed) *collapsed = cprev; 57 PetscFunctionReturn(0); 58 } 59 60 61 /* -------------------------------------------------------------------------- */ 62 /* 63 PCGAMGCreateGraph - create simple scaled scalar graph from matrix 64 65 Input Parameter: 66 . Amat - matrix 67 Output Parameter: 68 . a_Gmaat - eoutput scalar graph (symmetric?) 69 */ 70 #undef __FUNCT__ 71 #define __FUNCT__ "PCGAMGCreateGraph" 72 PetscErrorCode PCGAMGCreateGraph(Mat Amat, Mat *a_Gmat) 73 { 74 PetscErrorCode ierr; 75 PetscInt Istart,Iend,Ii,i,jj,kk,ncols,nloc,NN,MM,bs; 76 MPI_Comm comm; 77 Mat Gmat; 78 MatType mtype; 79 80 PetscFunctionBegin; 81 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 82 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 83 ierr = MatGetSize(Amat, &MM, &NN);CHKERRQ(ierr); 84 ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr); 85 nloc = (Iend-Istart)/bs; 86 87 #if defined PETSC_GAMG_USE_LOG 88 ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 89 #endif 90 91 if (bs > 1) { 92 const PetscScalar *vals; 93 const PetscInt *idx; 94 PetscInt *d_nnz, *o_nnz,*blockmask = NULL,maskcnt,*w0,*w1,*w2; 95 PetscBool ismpiaij,isseqaij; 96 97 /* 98 Determine the preallocation needed for the scalar matrix derived from the vector matrix. 99 */ 100 101 ierr = PetscObjectTypeCompare((PetscObject)Amat,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 102 ierr = PetscObjectTypeCompare((PetscObject)Amat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 103 104 ierr = PetscMalloc2(nloc, &d_nnz,isseqaij ? 0 : nloc, &o_nnz);CHKERRQ(ierr); 105 106 if (isseqaij) { 107 PetscInt max_d_nnz; 108 109 /* 110 Determine exact preallocation count for (sequential) scalar matrix 111 */ 112 ierr = MatSeqAIJGetMaxRowNonzeros(Amat,&max_d_nnz);CHKERRQ(ierr); 113 max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr); 114 ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr); 115 for (Ii = 0, jj = 0; Ii < Iend; Ii += bs, jj++) { 116 ierr = MatCollapseRows(Amat,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr); 117 } 118 ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr); 119 120 } else if (ismpiaij) { 121 Mat Daij,Oaij; 122 const PetscInt *garray; 123 PetscInt max_d_nnz; 124 125 ierr = MatMPIAIJGetSeqAIJ(Amat,&Daij,&Oaij,&garray);CHKERRQ(ierr); 126 127 /* 128 Determine exact preallocation count for diagonal block portion of scalar matrix 129 */ 130 ierr = MatSeqAIJGetMaxRowNonzeros(Daij,&max_d_nnz);CHKERRQ(ierr); 131 max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr); 132 ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr); 133 for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) { 134 ierr = MatCollapseRows(Daij,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr); 135 } 136 ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr); 137 138 /* 139 Over estimate (usually grossly over), preallocation count for off-diagonal portion of scalar matrix 140 */ 141 for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) { 142 o_nnz[jj] = 0; 143 for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */ 144 ierr = MatGetRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr); 145 o_nnz[jj] += ncols; 146 ierr = MatRestoreRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr); 147 } 148 if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc; 149 } 150 151 } else { 152 SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require AIJ matrix type"); 153 /* 154 155 This is O(nloc*nloc/bs) work! 156 157 This is accurate for the "diagonal" block of the matrix but will be grossly high for the 158 off diagonal block most of the time but could be too low for the off-diagonal. 159 160 This should be fixed to be accurate for the off-diagonal portion. Cannot just use a mask 161 for the off-diagonal portion since for huge matrices that would require too much memory per 162 MPI process. 163 */ 164 ierr = PetscMalloc1(nloc, &blockmask);CHKERRQ(ierr); 165 for (Ii = Istart, jj = 0; Ii < Iend; Ii += bs, jj++) { 166 o_nnz[jj] = 0; 167 ierr = PetscMemzero(blockmask,nloc*sizeof(PetscInt));CHKERRQ(ierr); 168 for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */ 169 ierr = MatGetRow(Amat,Ii+kk,&ncols,&idx,0);CHKERRQ(ierr); 170 for (i=0; i<ncols; i++) { 171 if (idx[i] >= Istart && idx[i] < Iend) { 172 blockmask[(idx[i] - Istart)/bs] = 1; 173 } 174 } 175 if (ncols > o_nnz[jj]) { 176 o_nnz[jj] = ncols; 177 if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc; 178 } 179 ierr = MatRestoreRow(Amat,Ii+kk,&ncols,&idx,0);CHKERRQ(ierr); 180 } 181 maskcnt = 0; 182 for (i=0; i<nloc; i++) { 183 if (blockmask[i]) maskcnt++; 184 } 185 d_nnz[jj] = maskcnt; 186 } 187 ierr = PetscFree(blockmask);CHKERRQ(ierr); 188 } 189 190 /* get scalar copy (norms) of matrix */ 191 ierr = MatGetType(Amat,&mtype);CHKERRQ(ierr); 192 ierr = MatCreate(comm, &Gmat);CHKERRQ(ierr); 193 ierr = MatSetSizes(Gmat,nloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 194 ierr = MatSetBlockSizes(Gmat, 1, 1);CHKERRQ(ierr); 195 ierr = MatSetType(Gmat, mtype);CHKERRQ(ierr); 196 ierr = MatSeqAIJSetPreallocation(Gmat,0,d_nnz);CHKERRQ(ierr); 197 ierr = MatMPIAIJSetPreallocation(Gmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 198 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 199 200 for (Ii = Istart; Ii < Iend; Ii++) { 201 PetscInt dest_row = Ii/bs; 202 ierr = MatGetRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 203 for (jj=0; jj<ncols; jj++) { 204 PetscInt dest_col = idx[jj]/bs; 205 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 206 ierr = MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);CHKERRQ(ierr); 207 } 208 ierr = MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 209 } 210 ierr = MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 211 ierr = MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 212 } else { 213 /* just copy scalar matrix - abs() not taken here but scaled later */ 214 ierr = MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);CHKERRQ(ierr); 215 } 216 217 #if defined PETSC_GAMG_USE_LOG 218 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 219 #endif 220 221 *a_Gmat = Gmat; 222 PetscFunctionReturn(0); 223 } 224 225 /* -------------------------------------------------------------------------- */ 226 /* 227 PCGAMGFilterGraph - filter (remove zero and possibly small values from the) graph and symetrize if needed 228 229 Input Parameter: 230 . vfilter - threshold paramter [0,1) 231 . symm - symetrize? 232 In/Output Parameter: 233 . a_Gmat - original graph 234 */ 235 #undef __FUNCT__ 236 #define __FUNCT__ "PCGAMGFilterGraph" 237 PetscErrorCode PCGAMGFilterGraph(Mat *a_Gmat,PetscReal vfilter,PetscBool symm) 238 { 239 PetscErrorCode ierr; 240 PetscInt Istart,Iend,Ii,jj,ncols,nnz0,nnz1, NN, MM, nloc; 241 PetscMPIInt rank; 242 Mat Gmat = *a_Gmat, tGmat, matTrans; 243 MPI_Comm comm; 244 const PetscScalar *vals; 245 const PetscInt *idx; 246 PetscInt *d_nnz, *o_nnz; 247 Vec diag; 248 MatType mtype; 249 250 PetscFunctionBegin; 251 #if defined PETSC_GAMG_USE_LOG 252 ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 253 #endif 254 /* scale Gmat for all values between -1 and 1 */ 255 ierr = MatCreateVecs(Gmat, &diag, 0);CHKERRQ(ierr); 256 ierr = MatGetDiagonal(Gmat, diag);CHKERRQ(ierr); 257 ierr = VecReciprocal(diag);CHKERRQ(ierr); 258 ierr = VecSqrtAbs(diag);CHKERRQ(ierr); 259 ierr = MatDiagonalScale(Gmat, diag, diag);CHKERRQ(ierr); 260 ierr = VecDestroy(&diag);CHKERRQ(ierr); 261 262 if (vfilter < 0.0 && !symm) { 263 /* Just use the provided matrix as the graph but make all values positive */ 264 MatInfo info; 265 PetscScalar *avals; 266 PetscBool isaij,ismpiaij; 267 ierr = PetscObjectTypeCompare((PetscObject)Gmat,MATSEQAIJ,&isaij);CHKERRQ(ierr); 268 ierr = PetscObjectTypeCompare((PetscObject)Gmat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 269 if (!isaij && !ismpiaij) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_USER,"Require (MPI)AIJ matrix type"); 270 if (isaij) { 271 ierr = MatGetInfo(Gmat,MAT_LOCAL,&info);CHKERRQ(ierr); 272 ierr = MatSeqAIJGetArray(Gmat,&avals);CHKERRQ(ierr); 273 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 274 ierr = MatSeqAIJRestoreArray(Gmat,&avals);CHKERRQ(ierr); 275 } else { 276 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)Gmat->data; 277 ierr = MatGetInfo(aij->A,MAT_LOCAL,&info);CHKERRQ(ierr); 278 ierr = MatSeqAIJGetArray(aij->A,&avals);CHKERRQ(ierr); 279 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 280 ierr = MatSeqAIJRestoreArray(aij->A,&avals);CHKERRQ(ierr); 281 ierr = MatGetInfo(aij->B,MAT_LOCAL,&info);CHKERRQ(ierr); 282 ierr = MatSeqAIJGetArray(aij->B,&avals);CHKERRQ(ierr); 283 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 284 ierr = MatSeqAIJRestoreArray(aij->B,&avals);CHKERRQ(ierr); 285 } 286 #if defined PETSC_GAMG_USE_LOG 287 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 288 #endif 289 PetscFunctionReturn(0); 290 } 291 292 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 293 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 294 ierr = MatGetOwnershipRange(Gmat, &Istart, &Iend);CHKERRQ(ierr); 295 nloc = Iend - Istart; 296 ierr = MatGetSize(Gmat, &MM, &NN);CHKERRQ(ierr); 297 298 if (symm) { 299 ierr = MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);CHKERRQ(ierr); 300 } 301 302 /* Determine upper bound on nonzeros needed in new filtered matrix */ 303 ierr = PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);CHKERRQ(ierr); 304 for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) { 305 ierr = MatGetRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 306 d_nnz[jj] = ncols; 307 o_nnz[jj] = ncols; 308 ierr = MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 309 if (symm) { 310 ierr = MatGetRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 311 d_nnz[jj] += ncols; 312 o_nnz[jj] += ncols; 313 ierr = MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 314 } 315 if (d_nnz[jj] > nloc) d_nnz[jj] = nloc; 316 if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc; 317 } 318 ierr = MatGetType(Gmat,&mtype);CHKERRQ(ierr); 319 ierr = MatCreate(comm, &tGmat);CHKERRQ(ierr); 320 ierr = MatSetSizes(tGmat,nloc,nloc,MM,MM);CHKERRQ(ierr); 321 ierr = MatSetBlockSizes(tGmat, 1, 1);CHKERRQ(ierr); 322 ierr = MatSetType(tGmat, mtype);CHKERRQ(ierr); 323 ierr = MatSeqAIJSetPreallocation(tGmat,0,d_nnz);CHKERRQ(ierr); 324 ierr = MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 325 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 326 if (symm) { 327 ierr = MatDestroy(&matTrans);CHKERRQ(ierr); 328 } else { 329 /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */ 330 ierr = MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 331 } 332 333 for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) { 334 ierr = MatGetRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 335 for (jj=0; jj<ncols; jj++,nnz0++) { 336 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 337 if (PetscRealPart(sv) > vfilter) { 338 nnz1++; 339 if (symm) { 340 sv *= 0.5; 341 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 342 ierr = MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);CHKERRQ(ierr); 343 } else { 344 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 345 } 346 } 347 } 348 ierr = MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 349 } 350 ierr = MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 351 ierr = MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 352 353 #if defined PETSC_GAMG_USE_LOG 354 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 355 #endif 356 357 #if defined(PETSC_USE_INFO) 358 { 359 double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc; 360 ierr = PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);CHKERRQ(ierr); 361 } 362 #endif 363 ierr = MatDestroy(&Gmat);CHKERRQ(ierr); 364 *a_Gmat = tGmat; 365 PetscFunctionReturn(0); 366 } 367 368 /* -------------------------------------------------------------------------- */ 369 /* 370 PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have npe > 1 371 372 Input Parameter: 373 . Gmat - MPIAIJ matrix for scattters 374 . data_sz - number of data terms per node (# cols in output) 375 . data_in[nloc*data_sz] - column oriented data 376 Output Parameter: 377 . a_stride - numbrt of rows of output 378 . a_data_out[stride*data_sz] - output data with ghosts 379 */ 380 #undef __FUNCT__ 381 #define __FUNCT__ "PCGAMGGetDataWithGhosts" 382 PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out) 383 { 384 PetscErrorCode ierr; 385 MPI_Comm comm; 386 Vec tmp_crds; 387 Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)Gmat->data; 388 PetscInt nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc; 389 PetscScalar *data_arr; 390 PetscReal *datas; 391 PetscBool isMPIAIJ; 392 393 PetscFunctionBegin; 394 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 395 ierr = PetscObjectTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);CHKERRQ(ierr); 396 ierr = MatGetOwnershipRange(Gmat, &my0, &Iend);CHKERRQ(ierr); 397 nloc = Iend - my0; 398 ierr = VecGetLocalSize(mpimat->lvec, &num_ghosts);CHKERRQ(ierr); 399 nnodes = num_ghosts + nloc; 400 *a_stride = nnodes; 401 ierr = MatCreateVecs(Gmat, &tmp_crds, 0);CHKERRQ(ierr); 402 403 ierr = PetscMalloc1(data_sz*nnodes, &datas);CHKERRQ(ierr); 404 for (dir=0; dir<data_sz; dir++) { 405 /* set local, and global */ 406 for (kk=0; kk<nloc; kk++) { 407 PetscInt gid = my0 + kk; 408 PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */ 409 datas[dir*nnodes + kk] = PetscRealPart(crd); 410 411 ierr = VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);CHKERRQ(ierr); 412 } 413 ierr = VecAssemblyBegin(tmp_crds);CHKERRQ(ierr); 414 ierr = VecAssemblyEnd(tmp_crds);CHKERRQ(ierr); 415 /* get ghost datas */ 416 ierr = VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 417 ierr = VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 418 ierr = VecGetArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 419 for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]); 420 ierr = VecRestoreArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 421 } 422 ierr = VecDestroy(&tmp_crds);CHKERRQ(ierr); 423 *a_data_out = datas; 424 PetscFunctionReturn(0); 425 } 426 427 428 /* 429 * 430 * PCGAMGHashTableCreate 431 */ 432 433 #undef __FUNCT__ 434 #define __FUNCT__ "PCGAMGHashTableCreate" 435 PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab) 436 { 437 PetscErrorCode ierr; 438 PetscInt kk; 439 440 PetscFunctionBegin; 441 a_tab->size = a_size; 442 ierr = PetscMalloc1(a_size, &a_tab->table);CHKERRQ(ierr); 443 ierr = PetscMalloc1(a_size, &a_tab->data);CHKERRQ(ierr); 444 for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1; 445 PetscFunctionReturn(0); 446 } 447 448 #undef __FUNCT__ 449 #define __FUNCT__ "PCGAMGHashTableDestroy" 450 PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab) 451 { 452 PetscErrorCode ierr; 453 454 PetscFunctionBegin; 455 ierr = PetscFree(a_tab->table);CHKERRQ(ierr); 456 ierr = PetscFree(a_tab->data);CHKERRQ(ierr); 457 PetscFunctionReturn(0); 458 } 459 460 #undef __FUNCT__ 461 #define __FUNCT__ "PCGAMGHashTableAdd" 462 PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data) 463 { 464 PetscInt kk,idx; 465 466 PetscFunctionBegin; 467 if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %d.",a_key); 468 for (kk = 0, idx = GAMG_HASH(a_key); 469 kk < a_tab->size; 470 kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) { 471 472 if (a_tab->table[idx] == a_key) { 473 /* exists */ 474 a_tab->data[idx] = a_data; 475 break; 476 } else if (a_tab->table[idx] == -1) { 477 /* add */ 478 a_tab->table[idx] = a_key; 479 a_tab->data[idx] = a_data; 480 break; 481 } 482 } 483 if (kk==a_tab->size) { 484 /* this is not to efficient, waiting until completely full */ 485 PetscInt oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data; 486 PetscErrorCode ierr; 487 488 a_tab->size = new_size; 489 490 ierr = PetscMalloc1(a_tab->size, &a_tab->table);CHKERRQ(ierr); 491 ierr = PetscMalloc1(a_tab->size, &a_tab->data);CHKERRQ(ierr); 492 493 for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1; 494 for (kk=0;kk<oldsize;kk++) { 495 if (oldtable[kk] != -1) { 496 ierr = PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]);CHKERRQ(ierr); 497 } 498 } 499 ierr = PetscFree(oldtable);CHKERRQ(ierr); 500 ierr = PetscFree(olddata);CHKERRQ(ierr); 501 ierr = PCGAMGHashTableAdd(a_tab, a_key, a_data);CHKERRQ(ierr); 502 } 503 PetscFunctionReturn(0); 504 } 505 506