1 2 #include <../src/mat/impls/aij/seq/aij.h> 3 4 /* 5 Redundant - remove later! 6 */ 7 #undef __FUNCT__ 8 #define __FUNCT__ "MatGetColumnIJ_SeqAIJ_Color_Redundant" 9 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color_Redundant(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 10 { 11 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 12 PetscErrorCode ierr; 13 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 14 PetscInt nz = a->i[m],row,*jj,mr,col; 15 PetscInt *cspidx; 16 17 PetscFunctionBegin; 18 *nn = n; 19 if (!ia) PetscFunctionReturn(0); 20 if (symmetric) { 21 SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"MatGetColumnIJ_SeqAIJ_Color() not supported for the case symmetric"); 22 ierr = MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,0,oshift,(PetscInt**)ia,(PetscInt**)ja);CHKERRQ(ierr); 23 } else { 24 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&collengths);CHKERRQ(ierr); 25 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 26 ierr = PetscMalloc((n+1)*sizeof(PetscInt),&cia);CHKERRQ(ierr); 27 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cja);CHKERRQ(ierr); 28 ierr = PetscMalloc((nz+1)*sizeof(PetscInt),&cspidx);CHKERRQ(ierr); 29 jj = a->j; 30 for (i=0; i<nz; i++) { 31 collengths[jj[i]]++; 32 } 33 cia[0] = oshift; 34 for (i=0; i<n; i++) { 35 cia[i+1] = cia[i] + collengths[i]; 36 } 37 ierr = PetscMemzero(collengths,n*sizeof(PetscInt));CHKERRQ(ierr); 38 jj = a->j; 39 for (row=0; row<m; row++) { 40 mr = a->i[row+1] - a->i[row]; 41 for (i=0; i<mr; i++) { 42 col = *jj++; 43 44 cspidx[cia[col] + collengths[col] - oshift] = a->i[row] + i; /* index of a->j */ 45 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 46 } 47 } 48 ierr = PetscFree(collengths);CHKERRQ(ierr); 49 *ia = cia; *ja = cja; 50 *spidx = cspidx; 51 } 52 PetscFunctionReturn(0); 53 } 54 55 #undef __FUNCT__ 56 #define __FUNCT__ "MatRestoreColumnIJ_SeqAIJ_Color_Redundant" 57 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color_Redundant(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 58 { 59 PetscErrorCode ierr; 60 61 PetscFunctionBegin; 62 if (!ia) PetscFunctionReturn(0); 63 64 ierr = PetscFree(*ia);CHKERRQ(ierr); 65 ierr = PetscFree(*ja);CHKERRQ(ierr); 66 ierr = PetscFree(*spidx);CHKERRQ(ierr); 67 PetscFunctionReturn(0); 68 } 69 70 /* -------------------------------------------*/ 71 72 #undef __FUNCT__ 73 #define __FUNCT__ "MatFDColoringCreate_SeqAIJ_den2sp" 74 /* 75 This routine optimizes MatFDColoringCreate_SeqAIJ() by using den2sp array 76 */ 77 PetscErrorCode MatFDColoringCreate_SeqAIJ_den2sp(Mat mat,ISColoring iscoloring,MatFDColoring c) 78 { 79 PetscErrorCode ierr; 80 PetscInt i,n,nrows,N,j,k,m,ncols,col; 81 const PetscInt *is,*rows,*ci,*cj; 82 PetscInt nis = iscoloring->n,*rowhit,*columnsforrow,l,bs = 1,*spidx,*spidxhit,*den2sp,*den2sp_i; 83 IS *isa; 84 PetscBool flg1,flg2; 85 Mat_SeqAIJ *csp = (Mat_SeqAIJ*)mat->data; 86 87 PetscFunctionBegin; 88 printf("MatFDColoringCreate_SeqAIJ_den2sp ...\n"); 89 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();"); 90 91 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 92 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 93 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 94 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 95 if (flg1 || flg2) { 96 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 97 } 98 99 N = mat->cmap->N/bs; 100 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 101 c->N = mat->cmap->N/bs; 102 c->m = mat->rmap->N/bs; 103 c->rstart = 0; 104 105 c->ncolors = nis; 106 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 107 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 108 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 109 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 110 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 111 ierr = PetscMalloc((csp->nz+1)*sizeof(PetscInt),&den2sp);CHKERRQ(ierr); 112 den2sp_i = den2sp; 113 114 ierr = MatGetColumnIJ_SeqAIJ_Color_Redundant(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 115 116 ierr = PetscMalloc((c->m+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 117 ierr = PetscMemzero(rowhit,c->m*sizeof(PetscInt));CHKERRQ(ierr); 118 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); // N=mat->cmap->N/bs or c->M ? 119 ierr = PetscMalloc((c->m+1)*sizeof(PetscInt),&spidxhit);CHKERRQ(ierr); 120 121 for (i=0; i<nis; i++) { /* loop over colors */ 122 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 123 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 124 125 c->ncolumns[i] = n; 126 if (n) { 127 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 128 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 129 } else { 130 c->columns[i] = 0; 131 } 132 133 /* fast, crude version requires O(N*N) work */ 134 for (j=0; j<n; j++) { /* loop over columns */ 135 col = is[j]; 136 rows = cj + ci[col]; 137 m = ci[col+1] - ci[col]; 138 for (k=0; k<m; k++) { /* loop over columns marking them in rowhit */ 139 rowhit[*rows] = col + 1; 140 spidxhit[*rows] = spidx[ci[col] + k]; 141 rows++; 142 } 143 } 144 /* get num of rows by counting the number of hits */ 145 nrows = 0; 146 for (j=0; j<N; j++) { 147 if (rowhit[j]) nrows++; 148 } 149 c->nrows[i] = nrows; 150 151 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 152 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 153 nrows = 0; 154 for (j=0; j<N; j++) { /* loop over rows */ 155 if (rowhit[j]) { 156 c->rows[i][nrows] = j; /* row index */ 157 c->columnsforrow[i][nrows] = rowhit[j] - 1; /* column index */ 158 den2sp_i[nrows++] = spidxhit[j]; 159 rowhit[j] = 0.0; /* zero rowhit for reuse */ 160 } 161 } 162 den2sp_i += nrows; 163 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 164 } 165 ierr = MatRestoreColumnIJ_SeqAIJ_Color_Redundant(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&spidx,NULL);CHKERRQ(ierr); 166 167 ierr = PetscFree(rowhit);CHKERRQ(ierr); 168 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 169 ierr = PetscFree(spidxhit);CHKERRQ(ierr); 170 171 /* Optimize by adding the vscale, and scaleforrow[][] fields -- needed for seqaij??? */ 172 /* 173 see the version for MPIAIJ 174 */ 175 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->rmap->n,PETSC_DETERMINE,0,NULL,&c->vscale);CHKERRQ(ierr); 176 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 177 for (k=0; k<c->ncolors; k++) { 178 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 179 for (l=0; l<c->nrows[k]; l++) { 180 col = c->columnsforrow[k][l]; 181 182 c->vscaleforrow[k][l] = col; 183 } 184 } 185 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 186 c->den2sp = den2sp; 187 mat->ops->fdcoloringapply = MatFDColoringApply_SeqAIJ; 188 PetscFunctionReturn(0); 189 } 190 191 /* --------------------------------------------------------------- */ 192 193 #undef __FUNCT__ 194 #define __FUNCT__ "MatFDColoringCreate_SeqAIJ" 195 /* 196 This routine is shared by AIJ and BAIJ matrices, since it operators only on the nonzero structure of the elements or blocks. 197 This is why it has the ugly code with the MatGetBlockSize() 198 */ 199 PetscErrorCode MatFDColoringCreate_SeqAIJ(Mat mat,ISColoring iscoloring,MatFDColoring c) 200 { 201 PetscErrorCode ierr; 202 PetscInt i,n,nrows,N,j,k,m,ncols,col; 203 const PetscInt *is,*rows,*ci,*cj; 204 PetscInt nis = iscoloring->n,*rowhit,*columnsforrow,l,bs = 1; 205 IS *isa; 206 PetscBool done,flg = PETSC_FALSE; 207 PetscBool flg1,flg2; 208 PetscBool new_impl=PETSC_FALSE; 209 210 PetscFunctionBegin; 211 printf("MatFDColoringCreate_SeqAIJ ...\n"); 212 ierr = PetscOptionsName("-new","using new impls","",&new_impl);CHKERRQ(ierr); 213 if (new_impl){ 214 ierr = MatFDColoringCreate_SeqAIJ_den2sp(mat,iscoloring,c);CHKERRQ(ierr); 215 PetscFunctionReturn(0); 216 } 217 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix must be assembled by calls to MatAssemblyBegin/End();"); 218 219 ierr = ISColoringGetIS(iscoloring,PETSC_IGNORE,&isa);CHKERRQ(ierr); 220 /* this is ugly way to get blocksize but cannot call MatGetBlockSize() because AIJ can have bs > 1 */ 221 ierr = PetscObjectTypeCompare((PetscObject)mat,MATSEQBAIJ,&flg1);CHKERRQ(ierr); 222 ierr = PetscObjectTypeCompare((PetscObject)mat,MATMPIBAIJ,&flg2);CHKERRQ(ierr); 223 if (flg1 || flg2) { 224 ierr = MatGetBlockSize(mat,&bs);CHKERRQ(ierr); 225 } 226 227 N = mat->cmap->N/bs; 228 c->M = mat->rmap->N/bs; /* set total rows, columns and local rows */ 229 c->N = mat->cmap->N/bs; 230 c->m = mat->rmap->N/bs; 231 c->rstart = 0; 232 233 c->ncolors = nis; 234 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->ncolumns);CHKERRQ(ierr); 235 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columns);CHKERRQ(ierr); 236 ierr = PetscMalloc(nis*sizeof(PetscInt),&c->nrows);CHKERRQ(ierr); 237 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->rows);CHKERRQ(ierr); 238 ierr = PetscMalloc(nis*sizeof(PetscInt*),&c->columnsforrow);CHKERRQ(ierr); 239 240 ierr = MatGetColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 241 if (!done) SETERRQ1(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"MatGetColumnIJ() not supported for matrix type %s",((PetscObject)mat)->type_name); 242 243 /* 244 Temporary option to allow for debugging/testing 245 */ 246 ierr = PetscOptionsGetBool(NULL,"-matfdcoloring_slow",&flg,NULL);CHKERRQ(ierr); 247 248 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&rowhit);CHKERRQ(ierr); 249 ierr = PetscMalloc((N+1)*sizeof(PetscInt),&columnsforrow);CHKERRQ(ierr); 250 251 for (i=0; i<nis; i++) { 252 ierr = ISGetLocalSize(isa[i],&n);CHKERRQ(ierr); 253 ierr = ISGetIndices(isa[i],&is);CHKERRQ(ierr); 254 255 c->ncolumns[i] = n; 256 if (n) { 257 ierr = PetscMalloc(n*sizeof(PetscInt),&c->columns[i]);CHKERRQ(ierr); 258 ierr = PetscMemcpy(c->columns[i],is,n*sizeof(PetscInt));CHKERRQ(ierr); 259 } else { 260 c->columns[i] = 0; 261 } 262 263 if (!flg) { /* ------------------------------------------------------------------------------*/ 264 /* fast, crude version requires O(N*N) work */ 265 ierr = PetscMemzero(rowhit,N*sizeof(PetscInt));CHKERRQ(ierr); 266 /* loop over columns*/ 267 for (j=0; j<n; j++) { 268 col = is[j]; 269 rows = cj + ci[col]; 270 m = ci[col+1] - ci[col]; 271 /* loop over columns marking them in rowhit */ 272 for (k=0; k<m; k++) { 273 rowhit[*rows++] = col + 1; 274 } 275 } 276 /* count the number of hits */ 277 nrows = 0; 278 for (j=0; j<N; j++) { 279 if (rowhit[j]) nrows++; 280 } 281 c->nrows[i] = nrows; 282 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 283 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 284 nrows = 0; 285 for (j=0; j<N; j++) { 286 if (rowhit[j]) { 287 c->rows[i][nrows] = j; 288 c->columnsforrow[i][nrows] = rowhit[j] - 1; 289 nrows++; 290 } 291 } 292 } else { /*-------------------------------------------------------------------------------*/ 293 /* slow version, using rowhit as a linked list */ 294 PetscInt currentcol,fm,mfm; 295 rowhit[N] = N; 296 nrows = 0; 297 /* loop over columns */ 298 for (j=0; j<n; j++) { 299 col = is[j]; 300 rows = cj + ci[col]; 301 m = ci[col+1] - ci[col]; 302 /* loop over columns marking them in rowhit */ 303 fm = N; /* fm points to first entry in linked list */ 304 for (k=0; k<m; k++) { 305 currentcol = *rows++; 306 /* is it already in the list? */ 307 do { 308 mfm = fm; 309 fm = rowhit[fm]; 310 } while (fm < currentcol); 311 /* not in list so add it */ 312 if (fm != currentcol) { 313 nrows++; 314 columnsforrow[currentcol] = col; 315 /* next three lines insert new entry into linked list */ 316 rowhit[mfm] = currentcol; 317 rowhit[currentcol] = fm; 318 fm = currentcol; 319 /* fm points to present position in list since we know the columns are sorted */ 320 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Detected invalid coloring"); 321 } 322 } 323 c->nrows[i] = nrows; 324 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->rows[i]);CHKERRQ(ierr); 325 ierr = PetscMalloc((nrows+1)*sizeof(PetscInt),&c->columnsforrow[i]);CHKERRQ(ierr); 326 /* now store the linked list of rows into c->rows[i] */ 327 nrows = 0; 328 fm = rowhit[N]; 329 do { 330 c->rows[i][nrows] = fm; 331 c->columnsforrow[i][nrows++] = columnsforrow[fm]; 332 fm = rowhit[fm]; 333 } while (fm < N); 334 } /* ---------------------------------------------------------------------------------------*/ 335 ierr = ISRestoreIndices(isa[i],&is);CHKERRQ(ierr); 336 } 337 ierr = MatRestoreColumnIJ(mat,0,PETSC_FALSE,PETSC_FALSE,&ncols,&ci,&cj,&done);CHKERRQ(ierr); 338 339 ierr = PetscFree(rowhit);CHKERRQ(ierr); 340 ierr = PetscFree(columnsforrow);CHKERRQ(ierr); 341 342 /* Optimize by adding the vscale, and scaleforrow[][] fields */ 343 /* 344 see the version for MPIAIJ 345 */ 346 ierr = VecCreateGhost(PetscObjectComm((PetscObject)mat),mat->rmap->n,PETSC_DETERMINE,0,NULL,&c->vscale);CHKERRQ(ierr); 347 ierr = PetscMalloc(c->ncolors*sizeof(PetscInt*),&c->vscaleforrow);CHKERRQ(ierr); 348 for (k=0; k<c->ncolors; k++) { 349 ierr = PetscMalloc((c->nrows[k]+1)*sizeof(PetscInt),&c->vscaleforrow[k]);CHKERRQ(ierr); 350 for (l=0; l<c->nrows[k]; l++) { 351 col = c->columnsforrow[k][l]; 352 353 c->vscaleforrow[k][l] = col; 354 } 355 } 356 ierr = ISColoringRestoreIS(iscoloring,&isa);CHKERRQ(ierr); 357 PetscFunctionReturn(0); 358 } 359 360