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 /* 153 154 This is O(nloc*nloc/bs) work! 155 156 This is accurate for the "diagonal" block of the matrix but will be grossly high for the 157 off diagonal block most of the time but could be too low for the off-diagonal. 158 159 This should be fixed to be accurate for the off-diagonal portion. Cannot just use a mask 160 for the off-diagonal portion since for huge matrices that would require too much memory per 161 MPI process. 162 */ 163 ierr = PetscMalloc1(nloc, &blockmask);CHKERRQ(ierr); 164 for (Ii = Istart, jj = 0; Ii < Iend; Ii += bs, jj++) { 165 o_nnz[jj] = 0; 166 ierr = PetscMemzero(blockmask,nloc*sizeof(PetscInt));CHKERRQ(ierr); 167 for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */ 168 ierr = MatGetRow(Amat,Ii+kk,&ncols,&idx,0);CHKERRQ(ierr); 169 for (i=0; i<ncols; i++) { 170 if (idx[i] >= Istart && idx[i] < Iend) { 171 blockmask[(idx[i] - Istart)/bs] = 1; 172 } 173 } 174 if (ncols > o_nnz[jj]) { 175 o_nnz[jj] = ncols; 176 if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc; 177 } 178 ierr = MatRestoreRow(Amat,Ii+kk,&ncols,&idx,0);CHKERRQ(ierr); 179 } 180 maskcnt = 0; 181 for (i=0; i<nloc; i++) { 182 if (blockmask[i]) maskcnt++; 183 } 184 d_nnz[jj] = maskcnt; 185 } 186 ierr = PetscFree(blockmask);CHKERRQ(ierr); 187 } 188 189 /* get scalar copy (norms) of matrix -- AIJ specific!!! */ 190 ierr = MatGetType(Amat,&mtype);CHKERRQ(ierr); 191 ierr = MatCreate(comm, &Gmat);CHKERRQ(ierr); 192 ierr = MatSetSizes(Gmat,nloc,nloc,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 193 ierr = MatSetBlockSizes(Gmat, 1, 1);CHKERRQ(ierr); 194 ierr = MatSetType(Gmat, mtype);CHKERRQ(ierr); 195 ierr = MatSeqAIJSetPreallocation(Gmat,0,d_nnz);CHKERRQ(ierr); 196 ierr = MatMPIAIJSetPreallocation(Gmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 197 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 198 199 for (Ii = Istart; Ii < Iend; Ii++) { 200 PetscInt dest_row = Ii/bs; 201 ierr = MatGetRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 202 for (jj=0; jj<ncols; jj++) { 203 PetscInt dest_col = idx[jj]/bs; 204 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 205 ierr = MatSetValues(Gmat,1,&dest_row,1,&dest_col,&sv,ADD_VALUES);CHKERRQ(ierr); 206 } 207 ierr = MatRestoreRow(Amat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 208 } 209 ierr = MatAssemblyBegin(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 210 ierr = MatAssemblyEnd(Gmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 211 } else { 212 /* just copy scalar matrix - abs() not taken here but scaled later */ 213 ierr = MatDuplicate(Amat, MAT_COPY_VALUES, &Gmat);CHKERRQ(ierr); 214 } 215 216 #if defined PETSC_GAMG_USE_LOG 217 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 218 #endif 219 220 *a_Gmat = Gmat; 221 PetscFunctionReturn(0); 222 } 223 224 /* -------------------------------------------------------------------------- */ 225 /* 226 PCGAMGFilterGraph - filter (remove zero and possibly small values from the) graph and symetrize if needed 227 228 Input Parameter: 229 . vfilter - threshold paramter [0,1) 230 . symm - symetrize? 231 In/Output Parameter: 232 . a_Gmat - original graph 233 */ 234 #undef __FUNCT__ 235 #define __FUNCT__ "PCGAMGFilterGraph" 236 PetscErrorCode PCGAMGFilterGraph(Mat *a_Gmat,PetscReal vfilter,PetscBool symm) 237 { 238 PetscErrorCode ierr; 239 PetscInt Istart,Iend,Ii,jj,ncols,nnz0,nnz1, NN, MM, nloc; 240 PetscMPIInt rank; 241 Mat Gmat = *a_Gmat, tGmat, matTrans; 242 MPI_Comm comm; 243 const PetscScalar *vals; 244 const PetscInt *idx; 245 PetscInt *d_nnz, *o_nnz; 246 Vec diag; 247 MatType mtype; 248 249 PetscFunctionBegin; 250 #if defined PETSC_GAMG_USE_LOG 251 ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 252 #endif 253 /* scale Gmat for all values between -1 and 1 */ 254 ierr = MatCreateVecs(Gmat, &diag, 0);CHKERRQ(ierr); 255 ierr = MatGetDiagonal(Gmat, diag);CHKERRQ(ierr); 256 ierr = VecReciprocal(diag);CHKERRQ(ierr); 257 ierr = VecSqrtAbs(diag);CHKERRQ(ierr); 258 ierr = MatDiagonalScale(Gmat, diag, diag);CHKERRQ(ierr); 259 ierr = VecDestroy(&diag);CHKERRQ(ierr); 260 261 if (vfilter < 0.0 && !symm) { 262 /* Just use the provided matrix as the graph but make all values positive */ 263 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)Gmat->data; 264 MatInfo info; 265 PetscScalar *avals; 266 PetscMPIInt size; 267 268 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)Gmat),&size);CHKERRQ(ierr); 269 if (size == 1) { 270 ierr = MatGetInfo(Gmat,MAT_LOCAL,&info);CHKERRQ(ierr); 271 ierr = MatSeqAIJGetArray(Gmat,&avals);CHKERRQ(ierr); 272 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 273 ierr = MatSeqAIJRestoreArray(Gmat,&avals);CHKERRQ(ierr); 274 } else { 275 ierr = MatGetInfo(aij->A,MAT_LOCAL,&info);CHKERRQ(ierr); 276 ierr = MatSeqAIJGetArray(aij->A,&avals);CHKERRQ(ierr); 277 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 278 ierr = MatSeqAIJRestoreArray(aij->A,&avals);CHKERRQ(ierr); 279 ierr = MatGetInfo(aij->B,MAT_LOCAL,&info);CHKERRQ(ierr); 280 ierr = MatSeqAIJGetArray(aij->B,&avals);CHKERRQ(ierr); 281 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 282 ierr = MatSeqAIJRestoreArray(aij->B,&avals);CHKERRQ(ierr); 283 } 284 #if defined PETSC_GAMG_USE_LOG 285 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 286 #endif 287 PetscFunctionReturn(0); 288 } 289 290 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 291 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 292 ierr = MatGetOwnershipRange(Gmat, &Istart, &Iend);CHKERRQ(ierr); 293 nloc = Iend - Istart; 294 ierr = MatGetSize(Gmat, &MM, &NN);CHKERRQ(ierr); 295 296 if (symm) { 297 ierr = MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);CHKERRQ(ierr); 298 } 299 300 /* Determine upper bound on nonzeros needed in new filtered matrix */ 301 ierr = PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);CHKERRQ(ierr); 302 for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) { 303 ierr = MatGetRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 304 d_nnz[jj] = ncols; 305 o_nnz[jj] = ncols; 306 ierr = MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 307 if (symm) { 308 ierr = MatGetRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 309 d_nnz[jj] += ncols; 310 o_nnz[jj] += ncols; 311 ierr = MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 312 } 313 if (d_nnz[jj] > nloc) d_nnz[jj] = nloc; 314 if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc; 315 } 316 ierr = MatGetType(Gmat,&mtype);CHKERRQ(ierr); 317 ierr = MatCreate(comm, &tGmat);CHKERRQ(ierr); 318 ierr = MatSetSizes(tGmat,nloc,nloc,MM,MM);CHKERRQ(ierr); 319 ierr = MatSetBlockSizes(tGmat, 1, 1);CHKERRQ(ierr); 320 ierr = MatSetType(tGmat, mtype);CHKERRQ(ierr); 321 ierr = MatSeqAIJSetPreallocation(tGmat,0,d_nnz);CHKERRQ(ierr); 322 ierr = MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 323 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 324 if (symm) { 325 ierr = MatDestroy(&matTrans);CHKERRQ(ierr); 326 } else { 327 /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */ 328 ierr = MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 329 } 330 331 for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) { 332 ierr = MatGetRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 333 for (jj=0; jj<ncols; jj++,nnz0++) { 334 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 335 if (PetscRealPart(sv) > vfilter) { 336 nnz1++; 337 if (symm) { 338 sv *= 0.5; 339 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 340 ierr = MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);CHKERRQ(ierr); 341 } else { 342 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 343 } 344 } 345 } 346 ierr = MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 347 } 348 ierr = MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 349 ierr = MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 350 351 #if defined PETSC_GAMG_USE_LOG 352 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 353 #endif 354 355 #if defined(PETSC_USE_INFO) 356 { 357 double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc; 358 ierr = PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);CHKERRQ(ierr); 359 } 360 #endif 361 ierr = MatDestroy(&Gmat);CHKERRQ(ierr); 362 *a_Gmat = tGmat; 363 PetscFunctionReturn(0); 364 } 365 366 /* -------------------------------------------------------------------------- */ 367 /* 368 PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have > 1 pe 369 370 Input Parameter: 371 . Gmat - MPIAIJ matrix for scattters 372 . data_sz - number of data terms per node (# cols in output) 373 . data_in[nloc*data_sz] - column oriented data 374 Output Parameter: 375 . a_stride - numbrt of rows of output 376 . a_data_out[stride*data_sz] - output data with ghosts 377 */ 378 #undef __FUNCT__ 379 #define __FUNCT__ "PCGAMGGetDataWithGhosts" 380 PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out) 381 { 382 PetscErrorCode ierr; 383 MPI_Comm comm; 384 Vec tmp_crds; 385 Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)Gmat->data; 386 PetscInt nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc; 387 PetscScalar *data_arr; 388 PetscReal *datas; 389 PetscBool isMPIAIJ; 390 391 PetscFunctionBegin; 392 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 393 ierr = PetscObjectTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);CHKERRQ(ierr); 394 ierr = MatGetOwnershipRange(Gmat, &my0, &Iend);CHKERRQ(ierr); 395 nloc = Iend - my0; 396 ierr = VecGetLocalSize(mpimat->lvec, &num_ghosts);CHKERRQ(ierr); 397 nnodes = num_ghosts + nloc; 398 *a_stride = nnodes; 399 ierr = MatCreateVecs(Gmat, &tmp_crds, 0);CHKERRQ(ierr); 400 401 ierr = PetscMalloc1(data_sz*nnodes, &datas);CHKERRQ(ierr); 402 for (dir=0; dir<data_sz; dir++) { 403 /* set local, and global */ 404 for (kk=0; kk<nloc; kk++) { 405 PetscInt gid = my0 + kk; 406 PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */ 407 datas[dir*nnodes + kk] = PetscRealPart(crd); 408 409 ierr = VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);CHKERRQ(ierr); 410 } 411 ierr = VecAssemblyBegin(tmp_crds);CHKERRQ(ierr); 412 ierr = VecAssemblyEnd(tmp_crds);CHKERRQ(ierr); 413 /* get ghost datas */ 414 ierr = VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 415 ierr = VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 416 ierr = VecGetArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 417 for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]); 418 ierr = VecRestoreArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 419 } 420 ierr = VecDestroy(&tmp_crds);CHKERRQ(ierr); 421 *a_data_out = datas; 422 PetscFunctionReturn(0); 423 } 424 425 426 /* 427 * 428 * PCGAMGHashTableCreate 429 */ 430 431 #undef __FUNCT__ 432 #define __FUNCT__ "PCGAMGHashTableCreate" 433 PetscErrorCode PCGAMGHashTableCreate(PetscInt a_size, PCGAMGHashTable *a_tab) 434 { 435 PetscErrorCode ierr; 436 PetscInt kk; 437 438 PetscFunctionBegin; 439 a_tab->size = a_size; 440 ierr = PetscMalloc1(a_size, &a_tab->table);CHKERRQ(ierr); 441 ierr = PetscMalloc1(a_size, &a_tab->data);CHKERRQ(ierr); 442 for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1; 443 PetscFunctionReturn(0); 444 } 445 446 #undef __FUNCT__ 447 #define __FUNCT__ "PCGAMGHashTableDestroy" 448 PetscErrorCode PCGAMGHashTableDestroy(PCGAMGHashTable *a_tab) 449 { 450 PetscErrorCode ierr; 451 452 PetscFunctionBegin; 453 ierr = PetscFree(a_tab->table);CHKERRQ(ierr); 454 ierr = PetscFree(a_tab->data);CHKERRQ(ierr); 455 PetscFunctionReturn(0); 456 } 457 458 #undef __FUNCT__ 459 #define __FUNCT__ "PCGAMGHashTableAdd" 460 PetscErrorCode PCGAMGHashTableAdd(PCGAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data) 461 { 462 PetscInt kk,idx; 463 464 PetscFunctionBegin; 465 if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %d.",a_key); 466 for (kk = 0, idx = GAMG_HASH(a_key); 467 kk < a_tab->size; 468 kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) { 469 470 if (a_tab->table[idx] == a_key) { 471 /* exists */ 472 a_tab->data[idx] = a_data; 473 break; 474 } else if (a_tab->table[idx] == -1) { 475 /* add */ 476 a_tab->table[idx] = a_key; 477 a_tab->data[idx] = a_data; 478 break; 479 } 480 } 481 if (kk==a_tab->size) { 482 /* this is not to efficient, waiting until completely full */ 483 PetscInt oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data; 484 PetscErrorCode ierr; 485 486 a_tab->size = new_size; 487 488 ierr = PetscMalloc1(a_tab->size, &a_tab->table);CHKERRQ(ierr); 489 ierr = PetscMalloc1(a_tab->size, &a_tab->data);CHKERRQ(ierr); 490 491 for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1; 492 for (kk=0;kk<oldsize;kk++) { 493 if (oldtable[kk] != -1) { 494 ierr = PCGAMGHashTableAdd(a_tab, oldtable[kk], olddata[kk]);CHKERRQ(ierr); 495 } 496 } 497 ierr = PetscFree(oldtable);CHKERRQ(ierr); 498 ierr = PetscFree(olddata);CHKERRQ(ierr); 499 ierr = PCGAMGHashTableAdd(a_tab, a_key, a_data);CHKERRQ(ierr); 500 } 501 PetscFunctionReturn(0); 502 } 503 504