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