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 #include <petsc/private/kspimpl.h> 7 8 #undef __FUNCT__ 9 #define __FUNCT__ "MatCollapseRow" 10 /* 11 Produces a set of block column indices of the matrix row, one for each block represented in the original row 12 13 n - the number of block indices in cc[] 14 cc - the block indices (must be large enough to contain the indices) 15 */ 16 PETSC_STATIC_INLINE PetscErrorCode MatCollapseRow(Mat Amat,PetscInt row,PetscInt bs,PetscInt *n,PetscInt *cc) 17 { 18 PetscInt cnt = -1,nidx,j; 19 const PetscInt *idx; 20 PetscErrorCode ierr; 21 22 PetscFunctionBegin; 23 ierr = MatGetRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr); 24 if (nidx) { 25 cnt = 0; 26 cc[cnt] = idx[0]/bs; 27 for (j=1; j<nidx; j++) { 28 if (cc[cnt] < idx[j]/bs) cc[++cnt] = idx[j]/bs; 29 } 30 } 31 ierr = MatRestoreRow(Amat,row,&nidx,&idx,NULL);CHKERRQ(ierr); 32 *n = cnt+1; 33 PetscFunctionReturn(0); 34 } 35 36 #undef __FUNCT__ 37 #define __FUNCT__ "MatCollapseRows" 38 /* 39 Produces a set of block column indices of the matrix block row, one for each block represented in the original set of rows 40 41 ncollapsed - the number of block indices 42 collapsed - the block indices (must be large enough to contain the indices) 43 */ 44 PETSC_STATIC_INLINE PetscErrorCode MatCollapseRows(Mat Amat,PetscInt start,PetscInt bs,PetscInt *w0,PetscInt *w1,PetscInt *w2,PetscInt *ncollapsed,PetscInt **collapsed) 45 { 46 PetscInt i,nprev,*cprev = w0,ncur = 0,*ccur = w1,*merged = w2,*cprevtmp; 47 PetscErrorCode ierr; 48 49 PetscFunctionBegin; 50 ierr = MatCollapseRow(Amat,start,bs,&nprev,cprev);CHKERRQ(ierr); 51 for (i=start+1; i<start+bs; i++) { 52 ierr = MatCollapseRow(Amat,i,bs,&ncur,ccur);CHKERRQ(ierr); 53 ierr = PetscMergeIntArray(nprev,cprev,ncur,ccur,&nprev,&merged);CHKERRQ(ierr); 54 cprevtmp = cprev; cprev = merged; merged = cprevtmp; 55 } 56 *ncollapsed = nprev; 57 if (collapsed) *collapsed = cprev; 58 PetscFunctionReturn(0); 59 } 60 61 62 /* -------------------------------------------------------------------------- */ 63 /* 64 PCGAMGCreateGraph - create simple scaled scalar graph from matrix 65 66 Input Parameter: 67 . Amat - matrix 68 Output Parameter: 69 . a_Gmaat - eoutput scalar graph (symmetric?) 70 */ 71 #undef __FUNCT__ 72 #define __FUNCT__ "PCGAMGCreateGraph" 73 PetscErrorCode PCGAMGCreateGraph(Mat Amat, Mat *a_Gmat) 74 { 75 PetscErrorCode ierr; 76 PetscInt Istart,Iend,Ii,i,jj,kk,ncols,nloc,NN,MM,bs; 77 MPI_Comm comm; 78 Mat Gmat; 79 MatType mtype; 80 81 PetscFunctionBegin; 82 ierr = PetscObjectGetComm((PetscObject)Amat,&comm);CHKERRQ(ierr); 83 ierr = MatGetOwnershipRange(Amat, &Istart, &Iend);CHKERRQ(ierr); 84 ierr = MatGetSize(Amat, &MM, &NN);CHKERRQ(ierr); 85 ierr = MatGetBlockSize(Amat, &bs);CHKERRQ(ierr); 86 nloc = (Iend-Istart)/bs; 87 88 #if defined PETSC_GAMG_USE_LOG 89 ierr = PetscLogEventBegin(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 90 #endif 91 92 if (bs > 1) { 93 const PetscScalar *vals; 94 const PetscInt *idx; 95 PetscInt *d_nnz, *o_nnz,*blockmask = NULL,maskcnt,*w0,*w1,*w2; 96 PetscBool ismpiaij,isseqaij; 97 98 /* 99 Determine the preallocation needed for the scalar matrix derived from the vector matrix. 100 */ 101 102 ierr = PetscObjectTypeCompare((PetscObject)Amat,MATSEQAIJ,&isseqaij);CHKERRQ(ierr); 103 ierr = PetscObjectTypeCompare((PetscObject)Amat,MATMPIAIJ,&ismpiaij);CHKERRQ(ierr); 104 105 ierr = PetscMalloc2(nloc, &d_nnz,isseqaij ? 0 : nloc, &o_nnz);CHKERRQ(ierr); 106 107 if (isseqaij) { 108 PetscInt max_d_nnz; 109 110 /* 111 Determine exact preallocation count for (sequential) scalar matrix 112 */ 113 ierr = MatSeqAIJGetMaxRowNonzeros(Amat,&max_d_nnz);CHKERRQ(ierr); 114 max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr); 115 ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr); 116 for (Ii = 0, jj = 0; Ii < Iend; Ii += bs, jj++) { 117 ierr = MatCollapseRows(Amat,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr); 118 } 119 ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr); 120 121 } else if (ismpiaij) { 122 Mat Daij,Oaij; 123 const PetscInt *garray; 124 PetscInt max_d_nnz; 125 126 ierr = MatMPIAIJGetSeqAIJ(Amat,&Daij,&Oaij,&garray);CHKERRQ(ierr); 127 128 /* 129 Determine exact preallocation count for diagonal block portion of scalar matrix 130 */ 131 ierr = MatSeqAIJGetMaxRowNonzeros(Daij,&max_d_nnz);CHKERRQ(ierr); 132 max_d_nnz = PetscMin(nloc,bs*max_d_nnz);CHKERRQ(ierr); 133 ierr = PetscMalloc3(max_d_nnz, &w0,max_d_nnz, &w1,max_d_nnz, &w2);CHKERRQ(ierr); 134 for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) { 135 ierr = MatCollapseRows(Daij,Ii,bs,w0,w1,w2,&d_nnz[jj],NULL);CHKERRQ(ierr); 136 } 137 ierr = PetscFree3(w0,w1,w2);CHKERRQ(ierr); 138 139 /* 140 Over estimate (usually grossly over), preallocation count for off-diagonal portion of scalar matrix 141 */ 142 for (Ii = 0, jj = 0; Ii < Iend - Istart; Ii += bs, jj++) { 143 o_nnz[jj] = 0; 144 for (kk=0; kk<bs; kk++) { /* rows that get collapsed to a single row */ 145 ierr = MatGetRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr); 146 o_nnz[jj] += ncols; 147 ierr = MatRestoreRow(Oaij,Ii+kk,&ncols,0,0);CHKERRQ(ierr); 148 } 149 if (o_nnz[jj] > (NN/bs-nloc)) o_nnz[jj] = NN/bs-nloc; 150 } 151 152 } else { 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 -- AIJ specific!!! */ 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 Mat_MPIAIJ *aij = (Mat_MPIAIJ*)Gmat->data; 265 MatInfo info; 266 PetscScalar *avals; 267 PetscMPIInt size; 268 269 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)Gmat),&size);CHKERRQ(ierr); 270 if (size == 1) { 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 ierr = MatGetInfo(aij->A,MAT_LOCAL,&info);CHKERRQ(ierr); 277 ierr = MatSeqAIJGetArray(aij->A,&avals);CHKERRQ(ierr); 278 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 279 ierr = MatSeqAIJRestoreArray(aij->A,&avals);CHKERRQ(ierr); 280 ierr = MatGetInfo(aij->B,MAT_LOCAL,&info);CHKERRQ(ierr); 281 ierr = MatSeqAIJGetArray(aij->B,&avals);CHKERRQ(ierr); 282 for (jj = 0; jj<info.nz_used; jj++) avals[jj] = PetscAbsScalar(avals[jj]); 283 ierr = MatSeqAIJRestoreArray(aij->B,&avals);CHKERRQ(ierr); 284 } 285 #if defined PETSC_GAMG_USE_LOG 286 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 287 #endif 288 PetscFunctionReturn(0); 289 } 290 291 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 292 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 293 ierr = MatGetOwnershipRange(Gmat, &Istart, &Iend);CHKERRQ(ierr); 294 nloc = Iend - Istart; 295 ierr = MatGetSize(Gmat, &MM, &NN);CHKERRQ(ierr); 296 297 if (symm) { 298 ierr = MatTranspose(Gmat, MAT_INITIAL_MATRIX, &matTrans);CHKERRQ(ierr); 299 } 300 301 /* Determine upper bound on nonzeros needed in new filtered matrix */ 302 ierr = PetscMalloc2(nloc, &d_nnz,nloc, &o_nnz);CHKERRQ(ierr); 303 for (Ii = Istart, jj = 0; Ii < Iend; Ii++, jj++) { 304 ierr = MatGetRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 305 d_nnz[jj] = ncols; 306 o_nnz[jj] = ncols; 307 ierr = MatRestoreRow(Gmat,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 308 if (symm) { 309 ierr = MatGetRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 310 d_nnz[jj] += ncols; 311 o_nnz[jj] += ncols; 312 ierr = MatRestoreRow(matTrans,Ii,&ncols,NULL,NULL);CHKERRQ(ierr); 313 } 314 if (d_nnz[jj] > nloc) d_nnz[jj] = nloc; 315 if (o_nnz[jj] > (MM-nloc)) o_nnz[jj] = MM - nloc; 316 } 317 ierr = MatGetType(Gmat,&mtype);CHKERRQ(ierr); 318 ierr = MatCreate(comm, &tGmat);CHKERRQ(ierr); 319 ierr = MatSetSizes(tGmat,nloc,nloc,MM,MM);CHKERRQ(ierr); 320 ierr = MatSetBlockSizes(tGmat, 1, 1);CHKERRQ(ierr); 321 ierr = MatSetType(tGmat, mtype);CHKERRQ(ierr); 322 ierr = MatSeqAIJSetPreallocation(tGmat,0,d_nnz);CHKERRQ(ierr); 323 ierr = MatMPIAIJSetPreallocation(tGmat,0,d_nnz,0,o_nnz);CHKERRQ(ierr); 324 ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr); 325 if (symm) { 326 ierr = MatDestroy(&matTrans);CHKERRQ(ierr); 327 } else { 328 /* all entries are generated locally so MatAssembly will be slightly faster for large process counts */ 329 ierr = MatSetOption(tGmat,MAT_NO_OFF_PROC_ENTRIES,PETSC_TRUE);CHKERRQ(ierr); 330 } 331 332 for (Ii = Istart, nnz0 = nnz1 = 0; Ii < Iend; Ii++) { 333 ierr = MatGetRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 334 for (jj=0; jj<ncols; jj++,nnz0++) { 335 PetscScalar sv = PetscAbs(PetscRealPart(vals[jj])); 336 if (PetscRealPart(sv) > vfilter) { 337 nnz1++; 338 if (symm) { 339 sv *= 0.5; 340 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 341 ierr = MatSetValues(tGmat,1,&idx[jj],1,&Ii,&sv,ADD_VALUES);CHKERRQ(ierr); 342 } else { 343 ierr = MatSetValues(tGmat,1,&Ii,1,&idx[jj],&sv,ADD_VALUES);CHKERRQ(ierr); 344 } 345 } 346 } 347 ierr = MatRestoreRow(Gmat,Ii,&ncols,&idx,&vals);CHKERRQ(ierr); 348 } 349 ierr = MatAssemblyBegin(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 350 ierr = MatAssemblyEnd(tGmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 351 352 #if defined PETSC_GAMG_USE_LOG 353 ierr = PetscLogEventEnd(petsc_gamg_setup_events[GRAPH],0,0,0,0);CHKERRQ(ierr); 354 #endif 355 356 #if defined(PETSC_USE_INFO) 357 { 358 double t1 = (!nnz0) ? 1. : 100.*(double)nnz1/(double)nnz0, t2 = (!nloc) ? 1. : (double)nnz0/(double)nloc; 359 ierr = PetscInfo4(*a_Gmat,"\t %g%% nnz after filtering, with threshold %g, %g nnz ave. (N=%D)\n",t1,vfilter,t2,MM);CHKERRQ(ierr); 360 } 361 #endif 362 ierr = MatDestroy(&Gmat);CHKERRQ(ierr); 363 *a_Gmat = tGmat; 364 PetscFunctionReturn(0); 365 } 366 367 /* -------------------------------------------------------------------------- */ 368 /* 369 PCGAMGGetDataWithGhosts - hacks into Mat MPIAIJ so this must have > 1 pe 370 371 Input Parameter: 372 . Gmat - MPIAIJ matrix for scattters 373 . data_sz - number of data terms per node (# cols in output) 374 . data_in[nloc*data_sz] - column oriented data 375 Output Parameter: 376 . a_stride - numbrt of rows of output 377 . a_data_out[stride*data_sz] - output data with ghosts 378 */ 379 #undef __FUNCT__ 380 #define __FUNCT__ "PCGAMGGetDataWithGhosts" 381 PetscErrorCode PCGAMGGetDataWithGhosts(Mat Gmat,PetscInt data_sz,PetscReal data_in[],PetscInt *a_stride,PetscReal **a_data_out) 382 { 383 PetscErrorCode ierr; 384 MPI_Comm comm; 385 Vec tmp_crds; 386 Mat_MPIAIJ *mpimat = (Mat_MPIAIJ*)Gmat->data; 387 PetscInt nnodes,num_ghosts,dir,kk,jj,my0,Iend,nloc; 388 PetscScalar *data_arr; 389 PetscReal *datas; 390 PetscBool isMPIAIJ; 391 392 PetscFunctionBegin; 393 ierr = PetscObjectGetComm((PetscObject)Gmat,&comm);CHKERRQ(ierr); 394 ierr = PetscObjectTypeCompare((PetscObject)Gmat, MATMPIAIJ, &isMPIAIJ);CHKERRQ(ierr); 395 ierr = MatGetOwnershipRange(Gmat, &my0, &Iend);CHKERRQ(ierr); 396 nloc = Iend - my0; 397 ierr = VecGetLocalSize(mpimat->lvec, &num_ghosts);CHKERRQ(ierr); 398 nnodes = num_ghosts + nloc; 399 *a_stride = nnodes; 400 ierr = MatCreateVecs(Gmat, &tmp_crds, 0);CHKERRQ(ierr); 401 402 ierr = PetscMalloc1(data_sz*nnodes, &datas);CHKERRQ(ierr); 403 for (dir=0; dir<data_sz; dir++) { 404 /* set local, and global */ 405 for (kk=0; kk<nloc; kk++) { 406 PetscInt gid = my0 + kk; 407 PetscScalar crd = (PetscScalar)data_in[dir*nloc + kk]; /* col oriented */ 408 datas[dir*nnodes + kk] = PetscRealPart(crd); 409 410 ierr = VecSetValues(tmp_crds, 1, &gid, &crd, INSERT_VALUES);CHKERRQ(ierr); 411 } 412 ierr = VecAssemblyBegin(tmp_crds);CHKERRQ(ierr); 413 ierr = VecAssemblyEnd(tmp_crds);CHKERRQ(ierr); 414 /* get ghost datas */ 415 ierr = VecScatterBegin(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 416 ierr = VecScatterEnd(mpimat->Mvctx,tmp_crds,mpimat->lvec,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 417 ierr = VecGetArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 418 for (kk=nloc,jj=0;jj<num_ghosts;kk++,jj++) datas[dir*nnodes + kk] = PetscRealPart(data_arr[jj]); 419 ierr = VecRestoreArray(mpimat->lvec, &data_arr);CHKERRQ(ierr); 420 } 421 ierr = VecDestroy(&tmp_crds);CHKERRQ(ierr); 422 *a_data_out = datas; 423 PetscFunctionReturn(0); 424 } 425 426 427 /* 428 * 429 * GAMGTableCreate 430 */ 431 432 #undef __FUNCT__ 433 #define __FUNCT__ "GAMGTableCreate" 434 PetscErrorCode GAMGTableCreate(PetscInt a_size, GAMGHashTable *a_tab) 435 { 436 PetscErrorCode ierr; 437 PetscInt kk; 438 439 PetscFunctionBegin; 440 a_tab->size = a_size; 441 ierr = PetscMalloc1(a_size, &a_tab->table);CHKERRQ(ierr); 442 ierr = PetscMalloc1(a_size, &a_tab->data);CHKERRQ(ierr); 443 for (kk=0; kk<a_size; kk++) a_tab->table[kk] = -1; 444 PetscFunctionReturn(0); 445 } 446 447 #undef __FUNCT__ 448 #define __FUNCT__ "GAMGTableDestroy" 449 PetscErrorCode GAMGTableDestroy(GAMGHashTable *a_tab) 450 { 451 PetscErrorCode ierr; 452 453 PetscFunctionBegin; 454 ierr = PetscFree(a_tab->table);CHKERRQ(ierr); 455 ierr = PetscFree(a_tab->data);CHKERRQ(ierr); 456 PetscFunctionReturn(0); 457 } 458 459 #undef __FUNCT__ 460 #define __FUNCT__ "GAMGTableAdd" 461 PetscErrorCode GAMGTableAdd(GAMGHashTable *a_tab, PetscInt a_key, PetscInt a_data) 462 { 463 PetscInt kk,idx; 464 465 PetscFunctionBegin; 466 if (a_key<0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_USER,"Negative key %d.",a_key); 467 for (kk = 0, idx = GAMG_HASH(a_key); 468 kk < a_tab->size; 469 kk++, idx = (idx==(a_tab->size-1)) ? 0 : idx + 1) { 470 471 if (a_tab->table[idx] == a_key) { 472 /* exists */ 473 a_tab->data[idx] = a_data; 474 break; 475 } else if (a_tab->table[idx] == -1) { 476 /* add */ 477 a_tab->table[idx] = a_key; 478 a_tab->data[idx] = a_data; 479 break; 480 } 481 } 482 if (kk==a_tab->size) { 483 /* this is not to efficient, waiting until completely full */ 484 PetscInt oldsize = a_tab->size, new_size = 2*a_tab->size + 5, *oldtable = a_tab->table, *olddata = a_tab->data; 485 PetscErrorCode ierr; 486 487 a_tab->size = new_size; 488 489 ierr = PetscMalloc1(a_tab->size, &a_tab->table);CHKERRQ(ierr); 490 ierr = PetscMalloc1(a_tab->size, &a_tab->data);CHKERRQ(ierr); 491 492 for (kk=0;kk<a_tab->size;kk++) a_tab->table[kk] = -1; 493 for (kk=0;kk<oldsize;kk++) { 494 if (oldtable[kk] != -1) { 495 ierr = GAMGTableAdd(a_tab, oldtable[kk], olddata[kk]);CHKERRQ(ierr); 496 } 497 } 498 ierr = PetscFree(oldtable);CHKERRQ(ierr); 499 ierr = PetscFree(olddata);CHKERRQ(ierr); 500 ierr = GAMGTableAdd(a_tab, a_key, a_data);CHKERRQ(ierr); 501 } 502 PetscFunctionReturn(0); 503 } 504 505