1 #define PETSCMAT_DLL 2 3 /* 4 Defines the basic matrix operations for the SBAIJ (compressed row) 5 matrix storage format. 6 */ 7 #include "../src/mat/impls/baij/seq/baij.h" /*I "petscmat.h" I*/ 8 #include "../src/mat/impls/sbaij/seq/sbaij.h" 9 10 extern PetscErrorCode MatSeqSBAIJSetNumericFactorization(Mat,PetscTruth); 11 #define CHUNKSIZE 10 12 13 /* 14 Checks for missing diagonals 15 */ 16 #undef __FUNCT__ 17 #define __FUNCT__ "MatMissingDiagonal_SeqSBAIJ" 18 PetscErrorCode MatMissingDiagonal_SeqSBAIJ(Mat A,PetscTruth *missing,PetscInt *dd) 19 { 20 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 21 PetscErrorCode ierr; 22 PetscInt *diag,*jj = a->j,i; 23 24 PetscFunctionBegin; 25 ierr = MatMarkDiagonal_SeqSBAIJ(A);CHKERRQ(ierr); 26 diag = a->diag; 27 *missing = PETSC_FALSE; 28 for (i=0; i<a->mbs; i++) { 29 if (jj[diag[i]] != i) { 30 *missing = PETSC_TRUE; 31 if (dd) *dd = i; 32 break; 33 } 34 } 35 PetscFunctionReturn(0); 36 } 37 38 #undef __FUNCT__ 39 #define __FUNCT__ "MatMarkDiagonal_SeqSBAIJ" 40 PetscErrorCode MatMarkDiagonal_SeqSBAIJ(Mat A) 41 { 42 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 43 PetscErrorCode ierr; 44 PetscInt i; 45 46 PetscFunctionBegin; 47 if (!a->diag) { 48 ierr = PetscMalloc(a->mbs*sizeof(PetscInt),&a->diag);CHKERRQ(ierr); 49 } 50 for (i=0; i<a->mbs; i++) a->diag[i] = a->i[i]; 51 PetscFunctionReturn(0); 52 } 53 54 #undef __FUNCT__ 55 #define __FUNCT__ "MatGetRowIJ_SeqSBAIJ" 56 static PetscErrorCode MatGetRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 57 { 58 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 59 PetscInt i,j,n = a->mbs,nz = a->i[n],bs = A->rmap->bs; 60 PetscErrorCode ierr; 61 62 PetscFunctionBegin; 63 *nn = n; 64 if (!ia) PetscFunctionReturn(0); 65 if (!blockcompressed) { 66 /* malloc & create the natural set of indices */ 67 ierr = PetscMalloc2((n+1)*bs,PetscInt,ia,nz*bs,PetscInt,ja);CHKERRQ(ierr); 68 for (i=0; i<n+1; i++) { 69 for (j=0; j<bs; j++) { 70 *ia[i*bs+j] = a->i[i]*bs+j+oshift; 71 } 72 } 73 for (i=0; i<nz; i++) { 74 for (j=0; j<bs; j++) { 75 *ja[i*bs+j] = a->j[i]*bs+j+oshift; 76 } 77 } 78 } else { /* blockcompressed */ 79 if (oshift == 1) { 80 /* temporarily add 1 to i and j indices */ 81 for (i=0; i<nz; i++) a->j[i]++; 82 for (i=0; i<n+1; i++) a->i[i]++; 83 } 84 *ia = a->i; *ja = a->j; 85 } 86 87 PetscFunctionReturn(0); 88 } 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "MatRestoreRowIJ_SeqSBAIJ" 92 static PetscErrorCode MatRestoreRowIJ_SeqSBAIJ(Mat A,PetscInt oshift,PetscTruth symmetric,PetscTruth blockcompressed,PetscInt *nn,PetscInt *ia[],PetscInt *ja[],PetscTruth *done) 93 { 94 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 95 PetscInt i,n = a->mbs,nz = a->i[n]; 96 PetscErrorCode ierr; 97 98 PetscFunctionBegin; 99 if (!ia) PetscFunctionReturn(0); 100 101 if (!blockcompressed) { 102 ierr = PetscFree2(*ia,*ja);CHKERRQ(ierr); 103 } else if (oshift == 1) { /* blockcompressed */ 104 for (i=0; i<nz; i++) a->j[i]--; 105 for (i=0; i<n+1; i++) a->i[i]--; 106 } 107 108 PetscFunctionReturn(0); 109 } 110 111 #undef __FUNCT__ 112 #define __FUNCT__ "MatDestroy_SeqSBAIJ" 113 PetscErrorCode MatDestroy_SeqSBAIJ(Mat A) 114 { 115 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 116 PetscErrorCode ierr; 117 118 PetscFunctionBegin; 119 #if defined(PETSC_USE_LOG) 120 PetscLogObjectState((PetscObject)A,"Rows=%D, NZ=%D",A->rmap->N,a->nz); 121 #endif 122 ierr = MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i);CHKERRQ(ierr); 123 if (a->row) {ierr = ISDestroy(a->row);CHKERRQ(ierr);} 124 if (a->col){ierr = ISDestroy(a->col);CHKERRQ(ierr);} 125 if (a->icol) {ierr = ISDestroy(a->icol);CHKERRQ(ierr);} 126 if (a->idiag) {ierr = PetscFree(a->idiag);CHKERRQ(ierr);} 127 if (a->inode.size) {ierr = PetscFree(a->inode.size);CHKERRQ(ierr);} 128 ierr = PetscFree(a->diag);CHKERRQ(ierr); 129 ierr = PetscFree2(a->imax,a->ilen);CHKERRQ(ierr); 130 ierr = PetscFree(a->solve_work);CHKERRQ(ierr); 131 ierr = PetscFree(a->relax_work);CHKERRQ(ierr); 132 ierr = PetscFree(a->solves_work);CHKERRQ(ierr); 133 ierr = PetscFree(a->mult_work);CHKERRQ(ierr); 134 ierr = PetscFree(a->saved_values);CHKERRQ(ierr); 135 ierr = PetscFree(a->xtoy);CHKERRQ(ierr); 136 137 ierr = PetscFree(a->inew);CHKERRQ(ierr); 138 ierr = PetscFree(a);CHKERRQ(ierr); 139 140 ierr = PetscObjectChangeTypeName((PetscObject)A,0);CHKERRQ(ierr); 141 ierr = PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C","",PETSC_NULL);CHKERRQ(ierr); 142 ierr = PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C","",PETSC_NULL);CHKERRQ(ierr); 143 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetColumnIndices_C","",PETSC_NULL);CHKERRQ(ierr); 144 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqaij_C","",PETSC_NULL);CHKERRQ(ierr); 145 ierr = PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqsbaij_seqbaij_C","",PETSC_NULL);CHKERRQ(ierr); 146 ierr = PetscObjectComposeFunction((PetscObject)A,"MatSeqSBAIJSetPreallocation_C","",PETSC_NULL);CHKERRQ(ierr); 147 PetscFunctionReturn(0); 148 } 149 150 #undef __FUNCT__ 151 #define __FUNCT__ "MatSetOption_SeqSBAIJ" 152 PetscErrorCode MatSetOption_SeqSBAIJ(Mat A,MatOption op,PetscTruth flg) 153 { 154 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 155 PetscErrorCode ierr; 156 157 PetscFunctionBegin; 158 switch (op) { 159 case MAT_ROW_ORIENTED: 160 a->roworiented = flg; 161 break; 162 case MAT_KEEP_NONZERO_PATTERN: 163 a->keepnonzeropattern = flg; 164 break; 165 case MAT_NEW_NONZERO_LOCATIONS: 166 a->nonew = (flg ? 0 : 1); 167 break; 168 case MAT_NEW_NONZERO_LOCATION_ERR: 169 a->nonew = (flg ? -1 : 0); 170 break; 171 case MAT_NEW_NONZERO_ALLOCATION_ERR: 172 a->nonew = (flg ? -2 : 0); 173 break; 174 case MAT_UNUSED_NONZERO_LOCATION_ERR: 175 a->nounused = (flg ? -1 : 0); 176 break; 177 case MAT_NEW_DIAGONALS: 178 case MAT_IGNORE_OFF_PROC_ENTRIES: 179 case MAT_USE_HASH_TABLE: 180 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 181 break; 182 case MAT_HERMITIAN: 183 if (flg) SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric"); 184 case MAT_SYMMETRIC: 185 case MAT_STRUCTURALLY_SYMMETRIC: 186 case MAT_SYMMETRY_ETERNAL: 187 if (!flg) SETERRQ(PETSC_ERR_SUP,"Matrix must be symmetric"); 188 ierr = PetscInfo1(A,"Option %s not relevent\n",MatOptions[op]);CHKERRQ(ierr); 189 break; 190 case MAT_IGNORE_LOWER_TRIANGULAR: 191 a->ignore_ltriangular = flg; 192 break; 193 case MAT_ERROR_LOWER_TRIANGULAR: 194 a->ignore_ltriangular = flg; 195 break; 196 case MAT_GETROW_UPPERTRIANGULAR: 197 a->getrow_utriangular = flg; 198 break; 199 default: 200 SETERRQ1(PETSC_ERR_SUP,"unknown option %d",op); 201 } 202 PetscFunctionReturn(0); 203 } 204 205 #undef __FUNCT__ 206 #define __FUNCT__ "MatGetRow_SeqSBAIJ" 207 PetscErrorCode MatGetRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *ncols,PetscInt **cols,PetscScalar **v) 208 { 209 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 210 PetscErrorCode ierr; 211 PetscInt itmp,i,j,k,M,*ai,*aj,bs,bn,bp,*cols_i,bs2; 212 MatScalar *aa,*aa_i; 213 PetscScalar *v_i; 214 215 PetscFunctionBegin; 216 if (A && !a->getrow_utriangular) SETERRQ(PETSC_ERR_SUP,"MatGetRow is not supported for SBAIJ matrix format. Getting the upper triangular part of row, run with -mat_getrow_uppertriangular, call MatSetOption(mat,MAT_GETROW_UPPERTRIANGULAR,PETSC_TRUE) or MatGetRowUpperTriangular()"); 217 /* Get the upper triangular part of the row */ 218 bs = A->rmap->bs; 219 ai = a->i; 220 aj = a->j; 221 aa = a->a; 222 bs2 = a->bs2; 223 224 if (row < 0 || row >= A->rmap->N) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE, "Row %D out of range", row); 225 226 bn = row/bs; /* Block number */ 227 bp = row % bs; /* Block position */ 228 M = ai[bn+1] - ai[bn]; 229 *ncols = bs*M; 230 231 if (v) { 232 *v = 0; 233 if (*ncols) { 234 ierr = PetscMalloc((*ncols+row)*sizeof(PetscScalar),v);CHKERRQ(ierr); 235 for (i=0; i<M; i++) { /* for each block in the block row */ 236 v_i = *v + i*bs; 237 aa_i = aa + bs2*(ai[bn] + i); 238 for (j=bp,k=0; j<bs2; j+=bs,k++) {v_i[k] = aa_i[j];} 239 } 240 } 241 } 242 243 if (cols) { 244 *cols = 0; 245 if (*ncols) { 246 ierr = PetscMalloc((*ncols+row)*sizeof(PetscInt),cols);CHKERRQ(ierr); 247 for (i=0; i<M; i++) { /* for each block in the block row */ 248 cols_i = *cols + i*bs; 249 itmp = bs*aj[ai[bn] + i]; 250 for (j=0; j<bs; j++) {cols_i[j] = itmp++;} 251 } 252 } 253 } 254 255 /*search column A(0:row-1,row) (=A(row,0:row-1)). Could be expensive! */ 256 /* this segment is currently removed, so only entries in the upper triangle are obtained */ 257 #ifdef column_search 258 v_i = *v + M*bs; 259 cols_i = *cols + M*bs; 260 for (i=0; i<bn; i++){ /* for each block row */ 261 M = ai[i+1] - ai[i]; 262 for (j=0; j<M; j++){ 263 itmp = aj[ai[i] + j]; /* block column value */ 264 if (itmp == bn){ 265 aa_i = aa + bs2*(ai[i] + j) + bs*bp; 266 for (k=0; k<bs; k++) { 267 *cols_i++ = i*bs+k; 268 *v_i++ = aa_i[k]; 269 } 270 *ncols += bs; 271 break; 272 } 273 } 274 } 275 #endif 276 PetscFunctionReturn(0); 277 } 278 279 #undef __FUNCT__ 280 #define __FUNCT__ "MatRestoreRow_SeqSBAIJ" 281 PetscErrorCode MatRestoreRow_SeqSBAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 282 { 283 PetscErrorCode ierr; 284 285 PetscFunctionBegin; 286 if (idx) {ierr = PetscFree(*idx);CHKERRQ(ierr);} 287 if (v) {ierr = PetscFree(*v);CHKERRQ(ierr);} 288 PetscFunctionReturn(0); 289 } 290 291 #undef __FUNCT__ 292 #define __FUNCT__ "MatGetRowUpperTriangular_SeqSBAIJ" 293 PetscErrorCode MatGetRowUpperTriangular_SeqSBAIJ(Mat A) 294 { 295 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 296 297 PetscFunctionBegin; 298 a->getrow_utriangular = PETSC_TRUE; 299 PetscFunctionReturn(0); 300 } 301 #undef __FUNCT__ 302 #define __FUNCT__ "MatRestoreRowUpperTriangular_SeqSBAIJ" 303 PetscErrorCode MatRestoreRowUpperTriangular_SeqSBAIJ(Mat A) 304 { 305 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 306 307 PetscFunctionBegin; 308 a->getrow_utriangular = PETSC_FALSE; 309 PetscFunctionReturn(0); 310 } 311 312 #undef __FUNCT__ 313 #define __FUNCT__ "MatTranspose_SeqSBAIJ" 314 PetscErrorCode MatTranspose_SeqSBAIJ(Mat A,MatReuse reuse,Mat *B) 315 { 316 PetscErrorCode ierr; 317 PetscFunctionBegin; 318 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 319 ierr = MatDuplicate(A,MAT_COPY_VALUES,B);CHKERRQ(ierr); 320 } 321 PetscFunctionReturn(0); 322 } 323 324 #undef __FUNCT__ 325 #define __FUNCT__ "MatView_SeqSBAIJ_ASCII" 326 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 327 { 328 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 329 PetscErrorCode ierr; 330 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 331 const char *name; 332 PetscViewerFormat format; 333 334 PetscFunctionBegin; 335 ierr = PetscObjectGetName((PetscObject)A,&name);CHKERRQ(ierr); 336 ierr = PetscViewerGetFormat(viewer,&format);CHKERRQ(ierr); 337 if (format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) { 338 ierr = PetscViewerASCIIPrintf(viewer," block size is %D\n",bs);CHKERRQ(ierr); 339 } else if (format == PETSC_VIEWER_ASCII_MATLAB) { 340 Mat aij; 341 342 if (A->factor && bs>1){ 343 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored with bs>1. MatView() with PETSC_VIEWER_ASCII_MATLAB is not supported and ignored!\n");CHKERRQ(ierr); 344 PetscFunctionReturn(0); 345 } 346 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&aij);CHKERRQ(ierr); 347 ierr = MatView(aij,viewer);CHKERRQ(ierr); 348 ierr = MatDestroy(aij);CHKERRQ(ierr); 349 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 350 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 351 for (i=0; i<a->mbs; i++) { 352 for (j=0; j<bs; j++) { 353 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 354 for (k=a->i[i]; k<a->i[i+1]; k++) { 355 for (l=0; l<bs; l++) { 356 #if defined(PETSC_USE_COMPLEX) 357 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 358 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 359 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 360 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0 && PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 361 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 362 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 363 } else if (PetscRealPart(a->a[bs2*k + l*bs + j]) != 0.0) { 364 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 365 } 366 #else 367 if (a->a[bs2*k + l*bs + j] != 0.0) { 368 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 369 } 370 #endif 371 } 372 } 373 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 374 } 375 } 376 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 377 } else if (format == PETSC_VIEWER_ASCII_FACTOR_INFO) { 378 PetscFunctionReturn(0); 379 } else { 380 if (A->factor && bs>1){ 381 ierr = PetscPrintf(PETSC_COMM_SELF,"Warning: matrix is factored. MatView_SeqSBAIJ_ASCII() may not display complete or logically correct entries!\n");CHKERRQ(ierr); 382 } 383 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_NO);CHKERRQ(ierr); 384 for (i=0; i<a->mbs; i++) { 385 for (j=0; j<bs; j++) { 386 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 387 for (k=a->i[i]; k<a->i[i+1]; k++) { 388 for (l=0; l<bs; l++) { 389 #if defined(PETSC_USE_COMPLEX) 390 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 391 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 392 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 393 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 394 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 395 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 396 } else { 397 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 398 } 399 #else 400 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 401 #endif 402 } 403 } 404 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 405 } 406 } 407 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_YES);CHKERRQ(ierr); 408 } 409 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 410 PetscFunctionReturn(0); 411 } 412 413 #undef __FUNCT__ 414 #define __FUNCT__ "MatView_SeqSBAIJ_Draw_Zoom" 415 static PetscErrorCode MatView_SeqSBAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 416 { 417 Mat A = (Mat) Aa; 418 Mat_SeqSBAIJ *a=(Mat_SeqSBAIJ*)A->data; 419 PetscErrorCode ierr; 420 PetscInt row,i,j,k,l,mbs=a->mbs,color,bs=A->rmap->bs,bs2=a->bs2; 421 PetscMPIInt rank; 422 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 423 MatScalar *aa; 424 MPI_Comm comm; 425 PetscViewer viewer; 426 427 PetscFunctionBegin; 428 /* 429 This is nasty. If this is called from an originally parallel matrix 430 then all processes call this,but only the first has the matrix so the 431 rest should return immediately. 432 */ 433 ierr = PetscObjectGetComm((PetscObject)draw,&comm);CHKERRQ(ierr); 434 ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr); 435 if (rank) PetscFunctionReturn(0); 436 437 ierr = PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer);CHKERRQ(ierr); 438 439 ierr = PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr);CHKERRQ(ierr); 440 PetscDrawString(draw, .3*(xl+xr), .3*(yl+yr), PETSC_DRAW_BLACK, "symmetric"); 441 442 /* loop over matrix elements drawing boxes */ 443 color = PETSC_DRAW_BLUE; 444 for (i=0,row=0; i<mbs; i++,row+=bs) { 445 for (j=a->i[i]; j<a->i[i+1]; j++) { 446 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 447 x_l = a->j[j]*bs; x_r = x_l + 1.0; 448 aa = a->a + j*bs2; 449 for (k=0; k<bs; k++) { 450 for (l=0; l<bs; l++) { 451 if (PetscRealPart(*aa++) >= 0.) continue; 452 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 453 } 454 } 455 } 456 } 457 color = PETSC_DRAW_CYAN; 458 for (i=0,row=0; i<mbs; i++,row+=bs) { 459 for (j=a->i[i]; j<a->i[i+1]; j++) { 460 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 461 x_l = a->j[j]*bs; x_r = x_l + 1.0; 462 aa = a->a + j*bs2; 463 for (k=0; k<bs; k++) { 464 for (l=0; l<bs; l++) { 465 if (PetscRealPart(*aa++) != 0.) continue; 466 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 467 } 468 } 469 } 470 } 471 472 color = PETSC_DRAW_RED; 473 for (i=0,row=0; i<mbs; i++,row+=bs) { 474 for (j=a->i[i]; j<a->i[i+1]; j++) { 475 y_l = A->rmap->N - row - 1.0; y_r = y_l + 1.0; 476 x_l = a->j[j]*bs; x_r = x_l + 1.0; 477 aa = a->a + j*bs2; 478 for (k=0; k<bs; k++) { 479 for (l=0; l<bs; l++) { 480 if (PetscRealPart(*aa++) <= 0.) continue; 481 ierr = PetscDrawRectangle(draw,x_l+k,y_l-l,x_r+k,y_r-l,color,color,color,color);CHKERRQ(ierr); 482 } 483 } 484 } 485 } 486 PetscFunctionReturn(0); 487 } 488 489 #undef __FUNCT__ 490 #define __FUNCT__ "MatView_SeqSBAIJ_Draw" 491 static PetscErrorCode MatView_SeqSBAIJ_Draw(Mat A,PetscViewer viewer) 492 { 493 PetscErrorCode ierr; 494 PetscReal xl,yl,xr,yr,w,h; 495 PetscDraw draw; 496 PetscTruth isnull; 497 498 PetscFunctionBegin; 499 ierr = PetscViewerDrawGetDraw(viewer,0,&draw);CHKERRQ(ierr); 500 ierr = PetscDrawIsNull(draw,&isnull);CHKERRQ(ierr); if (isnull) PetscFunctionReturn(0); 501 502 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer);CHKERRQ(ierr); 503 xr = A->rmap->N; yr = A->rmap->N; h = yr/10.0; w = xr/10.0; 504 xr += w; yr += h; xl = -w; yl = -h; 505 ierr = PetscDrawSetCoordinates(draw,xl,yl,xr,yr);CHKERRQ(ierr); 506 ierr = PetscDrawZoom(draw,MatView_SeqSBAIJ_Draw_Zoom,A);CHKERRQ(ierr); 507 ierr = PetscObjectCompose((PetscObject)A,"Zoomviewer",PETSC_NULL);CHKERRQ(ierr); 508 PetscFunctionReturn(0); 509 } 510 511 #undef __FUNCT__ 512 #define __FUNCT__ "MatView_SeqSBAIJ" 513 PetscErrorCode MatView_SeqSBAIJ(Mat A,PetscViewer viewer) 514 { 515 PetscErrorCode ierr; 516 PetscTruth iascii,isdraw; 517 518 PetscFunctionBegin; 519 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);CHKERRQ(ierr); 520 ierr = PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_DRAW,&isdraw);CHKERRQ(ierr); 521 if (iascii){ 522 ierr = MatView_SeqSBAIJ_ASCII(A,viewer);CHKERRQ(ierr); 523 } else if (isdraw) { 524 ierr = MatView_SeqSBAIJ_Draw(A,viewer);CHKERRQ(ierr); 525 } else { 526 Mat B; 527 ierr = MatConvert(A,MATSEQAIJ,MAT_INITIAL_MATRIX,&B);CHKERRQ(ierr); 528 ierr = MatView(B,viewer);CHKERRQ(ierr); 529 ierr = MatDestroy(B);CHKERRQ(ierr); 530 } 531 PetscFunctionReturn(0); 532 } 533 534 535 #undef __FUNCT__ 536 #define __FUNCT__ "MatGetValues_SeqSBAIJ" 537 PetscErrorCode MatGetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 538 { 539 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 540 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 541 PetscInt *ai = a->i,*ailen = a->ilen; 542 PetscInt brow,bcol,ridx,cidx,bs=A->rmap->bs,bs2=a->bs2; 543 MatScalar *ap,*aa = a->a; 544 545 PetscFunctionBegin; 546 for (k=0; k<m; k++) { /* loop over rows */ 547 row = im[k]; brow = row/bs; 548 if (row < 0) {v += n; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative row: %D",row); */ 549 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 550 rp = aj + ai[brow] ; ap = aa + bs2*ai[brow] ; 551 nrow = ailen[brow]; 552 for (l=0; l<n; l++) { /* loop over columns */ 553 if (in[l] < 0) {v++; continue;} /* SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Negative column: %D",in[l]); */ 554 if (in[l] >= A->cmap->n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->cmap->n-1); 555 col = in[l] ; 556 bcol = col/bs; 557 cidx = col%bs; 558 ridx = row%bs; 559 high = nrow; 560 low = 0; /* assume unsorted */ 561 while (high-low > 5) { 562 t = (low+high)/2; 563 if (rp[t] > bcol) high = t; 564 else low = t; 565 } 566 for (i=low; i<high; i++) { 567 if (rp[i] > bcol) break; 568 if (rp[i] == bcol) { 569 *v++ = ap[bs2*i+bs*cidx+ridx]; 570 goto finished; 571 } 572 } 573 *v++ = 0.0; 574 finished:; 575 } 576 } 577 PetscFunctionReturn(0); 578 } 579 580 581 #undef __FUNCT__ 582 #define __FUNCT__ "MatSetValuesBlocked_SeqSBAIJ" 583 PetscErrorCode MatSetValuesBlocked_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 584 { 585 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 586 PetscErrorCode ierr; 587 PetscInt *rp,k,low,high,t,ii,jj,row,nrow,i,col,l,rmax,N,lastcol = -1; 588 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen; 589 PetscInt *aj=a->j,nonew=a->nonew,bs2=a->bs2,bs=A->rmap->bs,stepval; 590 PetscTruth roworiented=a->roworiented; 591 const PetscScalar *value = v; 592 MatScalar *ap,*aa = a->a,*bap; 593 594 PetscFunctionBegin; 595 if (roworiented) { 596 stepval = (n-1)*bs; 597 } else { 598 stepval = (m-1)*bs; 599 } 600 for (k=0; k<m; k++) { /* loop over added rows */ 601 row = im[k]; 602 if (row < 0) continue; 603 #if defined(PETSC_USE_DEBUG) 604 if (row >= a->mbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 605 #endif 606 rp = aj + ai[row]; 607 ap = aa + bs2*ai[row]; 608 rmax = imax[row]; 609 nrow = ailen[row]; 610 low = 0; 611 high = nrow; 612 for (l=0; l<n; l++) { /* loop over added columns */ 613 if (in[l] < 0) continue; 614 col = in[l]; 615 #if defined(PETSC_USE_DEBUG) 616 if (col >= a->nbs) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 617 #endif 618 if (col < row) continue; /* ignore lower triangular block */ 619 if (roworiented) { 620 value = v + k*(stepval+bs)*bs + l*bs; 621 } else { 622 value = v + l*(stepval+bs)*bs + k*bs; 623 } 624 if (col <= lastcol) low = 0; else high = nrow; 625 lastcol = col; 626 while (high-low > 7) { 627 t = (low+high)/2; 628 if (rp[t] > col) high = t; 629 else low = t; 630 } 631 for (i=low; i<high; i++) { 632 if (rp[i] > col) break; 633 if (rp[i] == col) { 634 bap = ap + bs2*i; 635 if (roworiented) { 636 if (is == ADD_VALUES) { 637 for (ii=0; ii<bs; ii++,value+=stepval) { 638 for (jj=ii; jj<bs2; jj+=bs) { 639 bap[jj] += *value++; 640 } 641 } 642 } else { 643 for (ii=0; ii<bs; ii++,value+=stepval) { 644 for (jj=ii; jj<bs2; jj+=bs) { 645 bap[jj] = *value++; 646 } 647 } 648 } 649 } else { 650 if (is == ADD_VALUES) { 651 for (ii=0; ii<bs; ii++,value+=stepval) { 652 for (jj=0; jj<bs; jj++) { 653 *bap++ += *value++; 654 } 655 } 656 } else { 657 for (ii=0; ii<bs; ii++,value+=stepval) { 658 for (jj=0; jj<bs; jj++) { 659 *bap++ = *value++; 660 } 661 } 662 } 663 } 664 goto noinsert2; 665 } 666 } 667 if (nonew == 1) goto noinsert2; 668 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 669 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 670 N = nrow++ - 1; high++; 671 /* shift up all the later entries in this row */ 672 for (ii=N; ii>=i; ii--) { 673 rp[ii+1] = rp[ii]; 674 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 675 } 676 if (N >= i) { 677 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 678 } 679 rp[i] = col; 680 bap = ap + bs2*i; 681 if (roworiented) { 682 for (ii=0; ii<bs; ii++,value+=stepval) { 683 for (jj=ii; jj<bs2; jj+=bs) { 684 bap[jj] = *value++; 685 } 686 } 687 } else { 688 for (ii=0; ii<bs; ii++,value+=stepval) { 689 for (jj=0; jj<bs; jj++) { 690 *bap++ = *value++; 691 } 692 } 693 } 694 noinsert2:; 695 low = i; 696 } 697 ailen[row] = nrow; 698 } 699 PetscFunctionReturn(0); 700 } 701 702 #undef __FUNCT__ 703 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_Inode" 704 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_Inode(Mat A) 705 { 706 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 707 PetscErrorCode ierr; 708 const PetscInt *ai = a->i, *aj = a->j,*cols; 709 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 710 PetscTruth flag; 711 712 PetscFunctionBegin; 713 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 714 while (i < m){ 715 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 716 /* Limits the number of elements in a node to 'a->inode.limit' */ 717 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 718 nzy = ai[j+1] - ai[j]; 719 if (nzy != (nzx - j + i)) break; 720 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 721 if (!flag) break; 722 } 723 ns[node_count++] = blk_size; 724 i = j; 725 } 726 if (!a->inode.size && m && node_count > .9*m) { 727 ierr = PetscFree(ns);CHKERRQ(ierr); 728 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 729 } else { 730 a->inode.node_count = node_count; 731 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 732 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 733 ierr = PetscFree(ns);CHKERRQ(ierr); 734 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 735 736 /* count collections of adjacent columns in each inode */ 737 row = 0; 738 cnt = 0; 739 for (i=0; i<node_count; i++) { 740 cols = aj + ai[row] + a->inode.size[i]; 741 nz = ai[row+1] - ai[row] - a->inode.size[i]; 742 for (j=1; j<nz; j++) { 743 if (cols[j] != cols[j-1]+1) { 744 cnt++; 745 } 746 } 747 cnt++; 748 row += a->inode.size[i]; 749 } 750 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 751 cnt = 0; 752 row = 0; 753 for (i=0; i<node_count; i++) { 754 cols = aj + ai[row] + a->inode.size[i]; 755 CHKMEMQ; 756 counts[2*cnt] = cols[0]; 757 CHKMEMQ; 758 nz = ai[row+1] - ai[row] - a->inode.size[i]; 759 cnt2 = 1; 760 for (j=1; j<nz; j++) { 761 if (cols[j] != cols[j-1]+1) { 762 CHKMEMQ; 763 counts[2*(cnt++)+1] = cnt2; 764 counts[2*cnt] = cols[j]; 765 CHKMEMQ; 766 cnt2 = 1; 767 } else cnt2++; 768 } 769 CHKMEMQ; 770 counts[2*(cnt++)+1] = cnt2; 771 CHKMEMQ; 772 row += a->inode.size[i]; 773 } 774 ierr = PetscIntView(2*cnt,counts,0); 775 } 776 777 778 PetscFunctionReturn(0); 779 } 780 781 #undef __FUNCT__ 782 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 783 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 784 { 785 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 786 PetscErrorCode ierr; 787 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 788 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 789 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 790 MatScalar *aa = a->a,*ap; 791 792 PetscFunctionBegin; 793 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 794 795 if (m) rmax = ailen[0]; 796 for (i=1; i<mbs; i++) { 797 /* move each row back by the amount of empty slots (fshift) before it*/ 798 fshift += imax[i-1] - ailen[i-1]; 799 rmax = PetscMax(rmax,ailen[i]); 800 if (fshift) { 801 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 802 N = ailen[i]; 803 for (j=0; j<N; j++) { 804 ip[j-fshift] = ip[j]; 805 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 806 } 807 } 808 ai[i] = ai[i-1] + ailen[i-1]; 809 } 810 if (mbs) { 811 fshift += imax[mbs-1] - ailen[mbs-1]; 812 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 813 } 814 /* reset ilen and imax for each row */ 815 for (i=0; i<mbs; i++) { 816 ailen[i] = imax[i] = ai[i+1] - ai[i]; 817 } 818 a->nz = ai[mbs]; 819 820 /* diagonals may have moved, reset it */ 821 if (a->diag) { 822 ierr = PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 823 } 824 if (fshift && a->nounused == -1) { 825 SETERRQ4(PETSC_ERR_PLIB, "Unused space detected in matrix: %D X %D block size %D, %D unneeded", m, A->cmap->n, A->rmap->bs, fshift*bs2); 826 } 827 ierr = PetscInfo5(A,"Matrix size: %D X %D, block size %D; storage space: %D unneeded, %D used\n",m,A->rmap->N,A->rmap->bs,fshift*bs2,a->nz*bs2);CHKERRQ(ierr); 828 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 829 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 830 a->reallocs = 0; 831 A->info.nz_unneeded = (PetscReal)fshift*bs2; 832 a->idiagvalid = PETSC_FALSE; 833 if (bs2 == 1) { 834 ierr = MatAssemblyEnd_SeqSBAIJ_Inode(A);CHKERRQ(ierr); 835 } 836 PetscFunctionReturn(0); 837 } 838 839 /* 840 This function returns an array of flags which indicate the locations of contiguous 841 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 842 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 843 Assume: sizes should be long enough to hold all the values. 844 */ 845 #undef __FUNCT__ 846 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 847 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 848 { 849 PetscInt i,j,k,row; 850 PetscTruth flg; 851 852 PetscFunctionBegin; 853 for (i=0,j=0; i<n; j++) { 854 row = idx[i]; 855 if (row%bs!=0) { /* Not the begining of a block */ 856 sizes[j] = 1; 857 i++; 858 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 859 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 860 i++; 861 } else { /* Begining of the block, so check if the complete block exists */ 862 flg = PETSC_TRUE; 863 for (k=1; k<bs; k++) { 864 if (row+k != idx[i+k]) { /* break in the block */ 865 flg = PETSC_FALSE; 866 break; 867 } 868 } 869 if (flg) { /* No break in the bs */ 870 sizes[j] = bs; 871 i+= bs; 872 } else { 873 sizes[j] = 1; 874 i++; 875 } 876 } 877 } 878 *bs_max = j; 879 PetscFunctionReturn(0); 880 } 881 882 883 /* Only add/insert a(i,j) with i<=j (blocks). 884 Any a(i,j) with i>j input by user is ingored. 885 */ 886 887 #undef __FUNCT__ 888 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 889 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 890 { 891 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 892 PetscErrorCode ierr; 893 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 894 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 895 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 896 PetscInt ridx,cidx,bs2=a->bs2; 897 MatScalar *ap,value,*aa=a->a,*bap; 898 899 PetscFunctionBegin; 900 for (k=0; k<m; k++) { /* loop over added rows */ 901 row = im[k]; /* row number */ 902 brow = row/bs; /* block row number */ 903 if (row < 0) continue; 904 #if defined(PETSC_USE_DEBUG) 905 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 906 #endif 907 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 908 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 909 rmax = imax[brow]; /* maximum space allocated for this row */ 910 nrow = ailen[brow]; /* actual length of this row */ 911 low = 0; 912 913 for (l=0; l<n; l++) { /* loop over added columns */ 914 if (in[l] < 0) continue; 915 #if defined(PETSC_USE_DEBUG) 916 if (in[l] >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->N-1); 917 #endif 918 col = in[l]; 919 bcol = col/bs; /* block col number */ 920 921 if (brow > bcol) { 922 if (a->ignore_ltriangular){ 923 continue; /* ignore lower triangular values */ 924 } else { 925 SETERRQ(PETSC_ERR_USER,"Lower triangular value cannot be set for sbaij format. Ignoring these values, run with -mat_ignore_lower_triangular or call MatSetOption(mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_TRUE)"); 926 } 927 } 928 929 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 930 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 931 /* element value a(k,l) */ 932 if (roworiented) { 933 value = v[l + k*n]; 934 } else { 935 value = v[k + l*m]; 936 } 937 938 /* move pointer bap to a(k,l) quickly and add/insert value */ 939 if (col <= lastcol) low = 0; high = nrow; 940 lastcol = col; 941 while (high-low > 7) { 942 t = (low+high)/2; 943 if (rp[t] > bcol) high = t; 944 else low = t; 945 } 946 for (i=low; i<high; i++) { 947 if (rp[i] > bcol) break; 948 if (rp[i] == bcol) { 949 bap = ap + bs2*i + bs*cidx + ridx; 950 if (is == ADD_VALUES) *bap += value; 951 else *bap = value; 952 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 953 if (brow == bcol && ridx < cidx){ 954 bap = ap + bs2*i + bs*ridx + cidx; 955 if (is == ADD_VALUES) *bap += value; 956 else *bap = value; 957 } 958 goto noinsert1; 959 } 960 } 961 962 if (nonew == 1) goto noinsert1; 963 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 964 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 965 966 N = nrow++ - 1; high++; 967 /* shift up all the later entries in this row */ 968 for (ii=N; ii>=i; ii--) { 969 rp[ii+1] = rp[ii]; 970 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 971 } 972 if (N>=i) { 973 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 974 } 975 rp[i] = bcol; 976 ap[bs2*i + bs*cidx + ridx] = value; 977 noinsert1:; 978 low = i; 979 } 980 } /* end of loop over added columns */ 981 ailen[brow] = nrow; 982 } /* end of loop over added rows */ 983 PetscFunctionReturn(0); 984 } 985 986 #undef __FUNCT__ 987 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 988 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 989 { 990 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 991 Mat outA; 992 PetscErrorCode ierr; 993 PetscTruth row_identity; 994 995 PetscFunctionBegin; 996 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 997 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 998 if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported"); 999 if (inA->rmap->bs != 1) SETERRQ1(PETSC_ERR_SUP,"Matrix block size %D is not supported",inA->rmap->bs); /* Need to replace MatCholeskyFactorSymbolic_SeqSBAIJ_MSR()! */ 1000 1001 outA = inA; 1002 inA->factor = MAT_FACTOR_ICC; 1003 1004 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1005 ierr = MatSeqSBAIJSetNumericFactorization(inA,row_identity);CHKERRQ(ierr); 1006 1007 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1008 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr); } 1009 a->row = row; 1010 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1011 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr); } 1012 a->col = row; 1013 1014 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1015 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1016 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1017 1018 if (!a->solve_work) { 1019 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1020 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1021 } 1022 1023 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1024 PetscFunctionReturn(0); 1025 } 1026 1027 EXTERN_C_BEGIN 1028 #undef __FUNCT__ 1029 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1030 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1031 { 1032 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1033 PetscInt i,nz,n; 1034 1035 PetscFunctionBegin; 1036 nz = baij->maxnz; 1037 n = mat->cmap->n; 1038 for (i=0; i<nz; i++) { 1039 baij->j[i] = indices[i]; 1040 } 1041 baij->nz = nz; 1042 for (i=0; i<n; i++) { 1043 baij->ilen[i] = baij->imax[i]; 1044 } 1045 PetscFunctionReturn(0); 1046 } 1047 EXTERN_C_END 1048 1049 #undef __FUNCT__ 1050 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1051 /*@ 1052 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1053 in the matrix. 1054 1055 Input Parameters: 1056 + mat - the SeqSBAIJ matrix 1057 - indices - the column indices 1058 1059 Level: advanced 1060 1061 Notes: 1062 This can be called if you have precomputed the nonzero structure of the 1063 matrix and want to provide it to the matrix object to improve the performance 1064 of the MatSetValues() operation. 1065 1066 You MUST have set the correct numbers of nonzeros per row in the call to 1067 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1068 1069 MUST be called before any calls to MatSetValues() 1070 1071 .seealso: MatCreateSeqSBAIJ 1072 @*/ 1073 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1074 { 1075 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1076 1077 PetscFunctionBegin; 1078 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1079 PetscValidPointer(indices,2); 1080 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1081 if (f) { 1082 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1083 } else { 1084 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1085 } 1086 PetscFunctionReturn(0); 1087 } 1088 1089 #undef __FUNCT__ 1090 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1091 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1092 { 1093 PetscErrorCode ierr; 1094 1095 PetscFunctionBegin; 1096 /* If the two matrices have the same copy implementation, use fast copy. */ 1097 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1098 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1099 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1100 1101 if (a->i[A->rmap->N] != b->i[B->rmap->N]) { 1102 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1103 } 1104 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1105 } else { 1106 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1107 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1108 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1109 } 1110 PetscFunctionReturn(0); 1111 } 1112 1113 #undef __FUNCT__ 1114 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1115 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1116 { 1117 PetscErrorCode ierr; 1118 1119 PetscFunctionBegin; 1120 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1121 PetscFunctionReturn(0); 1122 } 1123 1124 #undef __FUNCT__ 1125 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1126 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1127 { 1128 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1129 PetscFunctionBegin; 1130 *array = a->a; 1131 PetscFunctionReturn(0); 1132 } 1133 1134 #undef __FUNCT__ 1135 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1136 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1137 { 1138 PetscFunctionBegin; 1139 PetscFunctionReturn(0); 1140 } 1141 1142 #include "petscblaslapack.h" 1143 #undef __FUNCT__ 1144 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1145 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1146 { 1147 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1148 PetscErrorCode ierr; 1149 PetscInt i,bs=Y->rmap->bs,bs2,j; 1150 PetscBLASInt one = 1,bnz = PetscBLASIntCast(x->nz); 1151 1152 PetscFunctionBegin; 1153 if (str == SAME_NONZERO_PATTERN) { 1154 PetscScalar alpha = a; 1155 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1156 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1157 if (y->xtoy && y->XtoY != X) { 1158 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1159 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1160 } 1161 if (!y->xtoy) { /* get xtoy */ 1162 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1163 y->XtoY = X; 1164 } 1165 bs2 = bs*bs; 1166 for (i=0; i<x->nz; i++) { 1167 j = 0; 1168 while (j < bs2){ 1169 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1170 j++; 1171 } 1172 } 1173 ierr = PetscInfo3(Y,"ratio of nnz_s(X)/nnz_s(Y): %D/%D = %G\n",bs2*x->nz,bs2*y->nz,(PetscReal)(bs2*x->nz)/(bs2*y->nz));CHKERRQ(ierr); 1174 } else { 1175 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1176 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1177 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1178 } 1179 PetscFunctionReturn(0); 1180 } 1181 1182 #undef __FUNCT__ 1183 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1184 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1185 { 1186 PetscFunctionBegin; 1187 *flg = PETSC_TRUE; 1188 PetscFunctionReturn(0); 1189 } 1190 1191 #undef __FUNCT__ 1192 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1193 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1194 { 1195 PetscFunctionBegin; 1196 *flg = PETSC_TRUE; 1197 PetscFunctionReturn(0); 1198 } 1199 1200 #undef __FUNCT__ 1201 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1202 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1203 { 1204 PetscFunctionBegin; 1205 *flg = PETSC_FALSE; 1206 PetscFunctionReturn(0); 1207 } 1208 1209 #undef __FUNCT__ 1210 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1211 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1212 { 1213 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1214 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1215 MatScalar *aa = a->a; 1216 1217 PetscFunctionBegin; 1218 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1219 PetscFunctionReturn(0); 1220 } 1221 1222 #undef __FUNCT__ 1223 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1224 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1225 { 1226 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1227 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1228 MatScalar *aa = a->a; 1229 1230 PetscFunctionBegin; 1231 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1232 PetscFunctionReturn(0); 1233 } 1234 1235 /* -------------------------------------------------------------------*/ 1236 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1237 MatGetRow_SeqSBAIJ, 1238 MatRestoreRow_SeqSBAIJ, 1239 MatMult_SeqSBAIJ_N, 1240 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1241 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1242 MatMultAdd_SeqSBAIJ_N, 1243 0, 1244 0, 1245 0, 1246 /*10*/ 0, 1247 0, 1248 MatCholeskyFactor_SeqSBAIJ, 1249 MatRelax_SeqSBAIJ, 1250 MatTranspose_SeqSBAIJ, 1251 /*15*/ MatGetInfo_SeqSBAIJ, 1252 MatEqual_SeqSBAIJ, 1253 MatGetDiagonal_SeqSBAIJ, 1254 MatDiagonalScale_SeqSBAIJ, 1255 MatNorm_SeqSBAIJ, 1256 /*20*/ 0, 1257 MatAssemblyEnd_SeqSBAIJ, 1258 MatSetOption_SeqSBAIJ, 1259 MatZeroEntries_SeqSBAIJ, 1260 /*24*/ 0, 1261 0, 1262 0, 1263 0, 1264 0, 1265 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1266 0, 1267 0, 1268 MatGetArray_SeqSBAIJ, 1269 MatRestoreArray_SeqSBAIJ, 1270 /*34*/ MatDuplicate_SeqSBAIJ, 1271 0, 1272 0, 1273 0, 1274 MatICCFactor_SeqSBAIJ, 1275 /*39*/ MatAXPY_SeqSBAIJ, 1276 MatGetSubMatrices_SeqSBAIJ, 1277 MatIncreaseOverlap_SeqSBAIJ, 1278 MatGetValues_SeqSBAIJ, 1279 MatCopy_SeqSBAIJ, 1280 /*44*/ 0, 1281 MatScale_SeqSBAIJ, 1282 0, 1283 0, 1284 0, 1285 /*49*/ 0, 1286 MatGetRowIJ_SeqSBAIJ, 1287 MatRestoreRowIJ_SeqSBAIJ, 1288 0, 1289 0, 1290 /*54*/ 0, 1291 0, 1292 0, 1293 0, 1294 MatSetValuesBlocked_SeqSBAIJ, 1295 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1296 0, 1297 0, 1298 0, 1299 0, 1300 /*64*/ 0, 1301 0, 1302 0, 1303 0, 1304 0, 1305 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1306 0, 1307 0, 1308 0, 1309 0, 1310 /*74*/ 0, 1311 0, 1312 0, 1313 0, 1314 0, 1315 /*79*/ 0, 1316 0, 1317 0, 1318 #if !defined(PETSC_USE_COMPLEX) 1319 MatGetInertia_SeqSBAIJ, 1320 #else 1321 0, 1322 #endif 1323 MatLoad_SeqSBAIJ, 1324 /*84*/ MatIsSymmetric_SeqSBAIJ, 1325 MatIsHermitian_SeqSBAIJ, 1326 MatIsStructurallySymmetric_SeqSBAIJ, 1327 0, 1328 0, 1329 /*89*/ 0, 1330 0, 1331 0, 1332 0, 1333 0, 1334 /*94*/ 0, 1335 0, 1336 0, 1337 0, 1338 0, 1339 /*99*/ 0, 1340 0, 1341 0, 1342 0, 1343 0, 1344 /*104*/0, 1345 MatRealPart_SeqSBAIJ, 1346 MatImaginaryPart_SeqSBAIJ, 1347 MatGetRowUpperTriangular_SeqSBAIJ, 1348 MatRestoreRowUpperTriangular_SeqSBAIJ, 1349 /*109*/0, 1350 0, 1351 0, 1352 0, 1353 MatMissingDiagonal_SeqSBAIJ 1354 }; 1355 1356 EXTERN_C_BEGIN 1357 #undef __FUNCT__ 1358 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1359 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1360 { 1361 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1362 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1363 PetscErrorCode ierr; 1364 1365 PetscFunctionBegin; 1366 if (aij->nonew != 1) { 1367 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1368 } 1369 1370 /* allocate space for values if not already there */ 1371 if (!aij->saved_values) { 1372 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1373 } 1374 1375 /* copy values over */ 1376 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1377 PetscFunctionReturn(0); 1378 } 1379 EXTERN_C_END 1380 1381 EXTERN_C_BEGIN 1382 #undef __FUNCT__ 1383 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1384 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1385 { 1386 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1387 PetscErrorCode ierr; 1388 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1389 1390 PetscFunctionBegin; 1391 if (aij->nonew != 1) { 1392 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1393 } 1394 if (!aij->saved_values) { 1395 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1396 } 1397 1398 /* copy values over */ 1399 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1400 PetscFunctionReturn(0); 1401 } 1402 EXTERN_C_END 1403 1404 EXTERN_C_BEGIN 1405 #undef __FUNCT__ 1406 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1407 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1408 { 1409 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1410 PetscErrorCode ierr; 1411 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1412 PetscTruth skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1413 1414 PetscFunctionBegin; 1415 B->preallocated = PETSC_TRUE; 1416 if (bs < 0) { 1417 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1418 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1419 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1420 bs = PetscAbs(bs); 1421 } 1422 if (nnz && newbs != bs) { 1423 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1424 } 1425 bs = newbs; 1426 1427 ierr = PetscMapSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1428 ierr = PetscMapSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1429 ierr = PetscMapSetUp(B->rmap);CHKERRQ(ierr); 1430 ierr = PetscMapSetUp(B->cmap);CHKERRQ(ierr); 1431 1432 mbs = B->rmap->N/bs; 1433 bs2 = bs*bs; 1434 1435 if (mbs*bs != B->rmap->N) { 1436 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1437 } 1438 1439 if (nz == MAT_SKIP_ALLOCATION) { 1440 skipallocation = PETSC_TRUE; 1441 nz = 0; 1442 } 1443 1444 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1445 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1446 if (nnz) { 1447 for (i=0; i<mbs; i++) { 1448 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1449 if (nnz[i] > mbs) SETERRQ3(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than block row length: local row %D value %D rowlength %D",i,nnz[i],mbs); 1450 } 1451 } 1452 1453 B->ops->mult = MatMult_SeqSBAIJ_N; 1454 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1455 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1456 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1457 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1458 if (!flg) { 1459 switch (bs) { 1460 case 1: 1461 B->ops->mult = MatMult_SeqSBAIJ_1; 1462 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1463 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1464 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1465 break; 1466 case 2: 1467 B->ops->mult = MatMult_SeqSBAIJ_2; 1468 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1469 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1470 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1471 break; 1472 case 3: 1473 B->ops->mult = MatMult_SeqSBAIJ_3; 1474 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1475 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1476 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1477 break; 1478 case 4: 1479 B->ops->mult = MatMult_SeqSBAIJ_4; 1480 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1481 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1482 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1483 break; 1484 case 5: 1485 B->ops->mult = MatMult_SeqSBAIJ_5; 1486 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1487 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1488 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1489 break; 1490 case 6: 1491 B->ops->mult = MatMult_SeqSBAIJ_6; 1492 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1493 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1494 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1495 break; 1496 case 7: 1497 B->ops->mult = MatMult_SeqSBAIJ_7; 1498 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1499 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1500 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1501 break; 1502 } 1503 } 1504 1505 b->mbs = mbs; 1506 b->nbs = mbs; 1507 if (!skipallocation) { 1508 if (!b->imax) { 1509 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1510 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1511 } 1512 if (!nnz) { 1513 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1514 else if (nz <= 0) nz = 1; 1515 for (i=0; i<mbs; i++) { 1516 b->imax[i] = nz; 1517 } 1518 nz = nz*mbs; /* total nz */ 1519 } else { 1520 nz = 0; 1521 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1522 } 1523 /* b->ilen will count nonzeros in each block row so far. */ 1524 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1525 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1526 1527 /* allocate the matrix space */ 1528 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1529 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1530 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1531 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1532 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1533 b->singlemalloc = PETSC_TRUE; 1534 1535 /* pointer to beginning of each row */ 1536 b->i[0] = 0; 1537 for (i=1; i<mbs+1; i++) { 1538 b->i[i] = b->i[i-1] + b->imax[i-1]; 1539 } 1540 b->free_a = PETSC_TRUE; 1541 b->free_ij = PETSC_TRUE; 1542 } else { 1543 b->free_a = PETSC_FALSE; 1544 b->free_ij = PETSC_FALSE; 1545 } 1546 1547 B->rmap->bs = bs; 1548 b->bs2 = bs2; 1549 b->nz = 0; 1550 b->maxnz = nz*bs2; 1551 1552 b->inew = 0; 1553 b->jnew = 0; 1554 b->anew = 0; 1555 b->a2anew = 0; 1556 b->permute = PETSC_FALSE; 1557 PetscFunctionReturn(0); 1558 } 1559 EXTERN_C_END 1560 1561 #undef __FUNCT__ 1562 #define __FUNCT__ "MatSeqBAIJSetNumericFactorization" 1563 /* 1564 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1565 */ 1566 PetscErrorCode MatSeqSBAIJSetNumericFactorization(Mat B,PetscTruth natural) 1567 { 1568 PetscErrorCode ierr; 1569 PetscTruth flg = PETSC_FALSE; 1570 PetscInt bs = B->rmap->bs; 1571 1572 PetscFunctionBegin; 1573 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1574 if (flg) bs = 8; 1575 1576 if (!natural) { 1577 switch (bs) { 1578 case 1: 1579 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1580 break; 1581 case 2: 1582 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1583 break; 1584 case 3: 1585 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1586 break; 1587 case 4: 1588 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1589 break; 1590 case 5: 1591 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1592 break; 1593 case 6: 1594 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1595 break; 1596 case 7: 1597 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1598 break; 1599 default: 1600 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1601 break; 1602 } 1603 } else { 1604 switch (bs) { 1605 case 1: 1606 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering; 1607 break; 1608 case 2: 1609 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1610 break; 1611 case 3: 1612 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1613 break; 1614 case 4: 1615 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1616 break; 1617 case 5: 1618 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1619 break; 1620 case 6: 1621 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1622 break; 1623 case 7: 1624 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1625 break; 1626 default: 1627 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1628 break; 1629 } 1630 } 1631 PetscFunctionReturn(0); 1632 } 1633 1634 EXTERN_C_BEGIN 1635 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1636 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1637 EXTERN_C_END 1638 1639 1640 EXTERN_C_BEGIN 1641 #undef __FUNCT__ 1642 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1643 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1644 { 1645 PetscInt n = A->rmap->n; 1646 PetscErrorCode ierr; 1647 1648 PetscFunctionBegin; 1649 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1650 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1651 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1652 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1653 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1654 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1655 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1656 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1657 (*B)->factor = ftype; 1658 PetscFunctionReturn(0); 1659 } 1660 EXTERN_C_END 1661 1662 EXTERN_C_BEGIN 1663 #undef __FUNCT__ 1664 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1665 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1666 { 1667 PetscFunctionBegin; 1668 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1669 *flg = PETSC_TRUE; 1670 } else { 1671 *flg = PETSC_FALSE; 1672 } 1673 PetscFunctionReturn(0); 1674 } 1675 EXTERN_C_END 1676 1677 EXTERN_C_BEGIN 1678 #if defined(PETSC_HAVE_MUMPS) 1679 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1680 #endif 1681 #if defined(PETSC_HAVE_SPOOLES) 1682 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1683 #endif 1684 #if defined(PETSC_HAVE_PASTIX) 1685 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1686 #endif 1687 EXTERN_C_END 1688 1689 /*MC 1690 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1691 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1692 1693 Options Database Keys: 1694 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1695 1696 Level: beginner 1697 1698 .seealso: MatCreateSeqSBAIJ 1699 M*/ 1700 1701 EXTERN_C_BEGIN 1702 #undef __FUNCT__ 1703 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1704 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1705 { 1706 Mat_SeqSBAIJ *b; 1707 PetscErrorCode ierr; 1708 PetscMPIInt size; 1709 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1710 1711 PetscFunctionBegin; 1712 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1713 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1714 1715 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1716 B->data = (void*)b; 1717 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1718 B->ops->destroy = MatDestroy_SeqSBAIJ; 1719 B->ops->view = MatView_SeqSBAIJ; 1720 B->mapping = 0; 1721 b->row = 0; 1722 b->icol = 0; 1723 b->reallocs = 0; 1724 b->saved_values = 0; 1725 b->inode.limit = 5; 1726 b->inode.max_limit = 5; 1727 1728 b->roworiented = PETSC_TRUE; 1729 b->nonew = 0; 1730 b->diag = 0; 1731 b->solve_work = 0; 1732 b->mult_work = 0; 1733 B->spptr = 0; 1734 b->keepnonzeropattern = PETSC_FALSE; 1735 b->xtoy = 0; 1736 b->XtoY = 0; 1737 1738 b->inew = 0; 1739 b->jnew = 0; 1740 b->anew = 0; 1741 b->a2anew = 0; 1742 b->permute = PETSC_FALSE; 1743 1744 b->ignore_ltriangular = PETSC_FALSE; 1745 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1746 1747 b->getrow_utriangular = PETSC_FALSE; 1748 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1749 1750 #if defined(PETSC_HAVE_PASTIX) 1751 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_pastix_C", 1752 "MatGetFactor_seqsbaij_pastix", 1753 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1754 #endif 1755 #if defined(PETSC_HAVE_SPOOLES) 1756 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_spooles_C", 1757 "MatGetFactor_seqsbaij_spooles", 1758 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1759 #endif 1760 #if defined(PETSC_HAVE_MUMPS) 1761 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_mumps_C", 1762 "MatGetFactor_seqsbaij_mumps", 1763 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1764 #endif 1765 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_seqsbaij_petsc_C", 1766 "MatGetFactorAvailable_seqsbaij_petsc", 1767 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1768 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_petsc_C", 1769 "MatGetFactor_seqsbaij_petsc", 1770 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1771 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1772 "MatStoreValues_SeqSBAIJ", 1773 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1775 "MatRetrieveValues_SeqSBAIJ", 1776 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1777 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1778 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1779 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1781 "MatConvert_SeqSBAIJ_SeqAIJ", 1782 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1783 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1784 "MatConvert_SeqSBAIJ_SeqBAIJ", 1785 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1786 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1787 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1788 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1789 1790 B->symmetric = PETSC_TRUE; 1791 B->structurally_symmetric = PETSC_TRUE; 1792 B->symmetric_set = PETSC_TRUE; 1793 B->structurally_symmetric_set = PETSC_TRUE; 1794 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1795 1796 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1797 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1798 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1799 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1800 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1801 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",PETSC_NULL,b->inode.limit,&b->inode.limit,PETSC_NULL);CHKERRQ(ierr); 1802 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1803 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1804 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1805 1806 PetscFunctionReturn(0); 1807 } 1808 EXTERN_C_END 1809 1810 #undef __FUNCT__ 1811 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1812 /*@C 1813 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1814 compressed row) format. For good matrix assembly performance the 1815 user should preallocate the matrix storage by setting the parameter nz 1816 (or the array nnz). By setting these parameters accurately, performance 1817 during matrix assembly can be increased by more than a factor of 50. 1818 1819 Collective on Mat 1820 1821 Input Parameters: 1822 + A - the symmetric matrix 1823 . bs - size of block 1824 . nz - number of block nonzeros per block row (same for all rows) 1825 - nnz - array containing the number of block nonzeros in the upper triangular plus 1826 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1827 1828 Options Database Keys: 1829 . -mat_no_unroll - uses code that does not unroll the loops in the 1830 block calculations (much slower) 1831 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1832 1833 Level: intermediate 1834 1835 Notes: 1836 Specify the preallocated storage with either nz or nnz (not both). 1837 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1838 allocation. For additional details, see the users manual chapter on 1839 matrices. 1840 1841 You can call MatGetInfo() to get information on how effective the preallocation was; 1842 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1843 You can also run with the option -info and look for messages with the string 1844 malloc in them to see if additional memory allocation was needed. 1845 1846 If the nnz parameter is given then the nz parameter is ignored 1847 1848 1849 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1850 @*/ 1851 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1852 { 1853 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1854 1855 PetscFunctionBegin; 1856 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1857 if (f) { 1858 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1859 } 1860 PetscFunctionReturn(0); 1861 } 1862 1863 #undef __FUNCT__ 1864 #define __FUNCT__ "MatCreateSeqSBAIJ" 1865 /*@C 1866 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1867 compressed row) format. For good matrix assembly performance the 1868 user should preallocate the matrix storage by setting the parameter nz 1869 (or the array nnz). By setting these parameters accurately, performance 1870 during matrix assembly can be increased by more than a factor of 50. 1871 1872 Collective on MPI_Comm 1873 1874 Input Parameters: 1875 + comm - MPI communicator, set to PETSC_COMM_SELF 1876 . bs - size of block 1877 . m - number of rows, or number of columns 1878 . nz - number of block nonzeros per block row (same for all rows) 1879 - nnz - array containing the number of block nonzeros in the upper triangular plus 1880 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1881 1882 Output Parameter: 1883 . A - the symmetric matrix 1884 1885 Options Database Keys: 1886 . -mat_no_unroll - uses code that does not unroll the loops in the 1887 block calculations (much slower) 1888 . -mat_block_size - size of the blocks to use 1889 1890 Level: intermediate 1891 1892 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1893 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1894 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1895 1896 Notes: 1897 The number of rows and columns must be divisible by blocksize. 1898 This matrix type does not support complex Hermitian operation. 1899 1900 Specify the preallocated storage with either nz or nnz (not both). 1901 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1902 allocation. For additional details, see the users manual chapter on 1903 matrices. 1904 1905 If the nnz parameter is given then the nz parameter is ignored 1906 1907 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1908 @*/ 1909 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1910 { 1911 PetscErrorCode ierr; 1912 1913 PetscFunctionBegin; 1914 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1915 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1916 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1917 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1918 PetscFunctionReturn(0); 1919 } 1920 1921 #undef __FUNCT__ 1922 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1923 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1924 { 1925 Mat C; 1926 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1927 PetscErrorCode ierr; 1928 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1929 1930 PetscFunctionBegin; 1931 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1932 1933 *B = 0; 1934 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1935 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1936 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1937 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1938 c = (Mat_SeqSBAIJ*)C->data; 1939 1940 C->preallocated = PETSC_TRUE; 1941 C->factor = A->factor; 1942 c->row = 0; 1943 c->icol = 0; 1944 c->saved_values = 0; 1945 c->keepnonzeropattern = a->keepnonzeropattern; 1946 C->assembled = PETSC_TRUE; 1947 1948 ierr = PetscMapCopy(((PetscObject)A)->comm,A->rmap,C->rmap);CHKERRQ(ierr); 1949 ierr = PetscMapCopy(((PetscObject)A)->comm,A->cmap,C->cmap);CHKERRQ(ierr); 1950 c->bs2 = a->bs2; 1951 c->mbs = a->mbs; 1952 c->nbs = a->nbs; 1953 1954 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1955 for (i=0; i<mbs; i++) { 1956 c->imax[i] = a->imax[i]; 1957 c->ilen[i] = a->ilen[i]; 1958 } 1959 1960 /* allocate the matrix space */ 1961 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1962 c->singlemalloc = PETSC_TRUE; 1963 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1964 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1965 if (mbs > 0) { 1966 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1967 if (cpvalues == MAT_COPY_VALUES) { 1968 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1969 } else { 1970 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1971 } 1972 } 1973 1974 c->roworiented = a->roworiented; 1975 c->nonew = a->nonew; 1976 1977 if (a->diag) { 1978 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1979 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1980 for (i=0; i<mbs; i++) { 1981 c->diag[i] = a->diag[i]; 1982 } 1983 } else c->diag = 0; 1984 c->nz = a->nz; 1985 c->maxnz = a->maxnz; 1986 c->solve_work = 0; 1987 c->mult_work = 0; 1988 c->free_a = PETSC_TRUE; 1989 c->free_ij = PETSC_TRUE; 1990 *B = C; 1991 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 1992 PetscFunctionReturn(0); 1993 } 1994 1995 #undef __FUNCT__ 1996 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1997 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 1998 { 1999 Mat_SeqSBAIJ *a; 2000 Mat B; 2001 PetscErrorCode ierr; 2002 int fd; 2003 PetscMPIInt size; 2004 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2005 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2006 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2007 PetscInt *masked,nmask,tmp,bs2,ishift; 2008 PetscScalar *aa; 2009 MPI_Comm comm = ((PetscObject)viewer)->comm; 2010 2011 PetscFunctionBegin; 2012 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2013 bs2 = bs*bs; 2014 2015 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2016 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 2017 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2018 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2019 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2020 M = header[1]; N = header[2]; nz = header[3]; 2021 2022 if (header[3] < 0) { 2023 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2024 } 2025 2026 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 2027 2028 /* 2029 This code adds extra rows to make sure the number of rows is 2030 divisible by the blocksize 2031 */ 2032 mbs = M/bs; 2033 extra_rows = bs - M + bs*(mbs); 2034 if (extra_rows == bs) extra_rows = 0; 2035 else mbs++; 2036 if (extra_rows) { 2037 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2038 } 2039 2040 /* read in row lengths */ 2041 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2042 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2043 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2044 2045 /* read in column indices */ 2046 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2047 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2048 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2049 2050 /* loop over row lengths determining block row lengths */ 2051 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2052 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2053 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 2054 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2055 masked = mask + mbs; 2056 rowcount = 0; nzcount = 0; 2057 for (i=0; i<mbs; i++) { 2058 nmask = 0; 2059 for (j=0; j<bs; j++) { 2060 kmax = rowlengths[rowcount]; 2061 for (k=0; k<kmax; k++) { 2062 tmp = jj[nzcount++]/bs; /* block col. index */ 2063 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2064 } 2065 rowcount++; 2066 } 2067 s_browlengths[i] += nmask; 2068 2069 /* zero out the mask elements we set */ 2070 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2071 } 2072 2073 /* create our matrix */ 2074 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2075 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2076 ierr = MatSetType(B,type);CHKERRQ(ierr); 2077 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2078 a = (Mat_SeqSBAIJ*)B->data; 2079 2080 /* set matrix "i" values */ 2081 a->i[0] = 0; 2082 for (i=1; i<= mbs; i++) { 2083 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2084 a->ilen[i-1] = s_browlengths[i-1]; 2085 } 2086 a->nz = a->i[mbs]; 2087 2088 /* read in nonzero values */ 2089 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2090 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2091 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2092 2093 /* set "a" and "j" values into matrix */ 2094 nzcount = 0; jcount = 0; 2095 for (i=0; i<mbs; i++) { 2096 nzcountb = nzcount; 2097 nmask = 0; 2098 for (j=0; j<bs; j++) { 2099 kmax = rowlengths[i*bs+j]; 2100 for (k=0; k<kmax; k++) { 2101 tmp = jj[nzcount++]/bs; /* block col. index */ 2102 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2103 } 2104 } 2105 /* sort the masked values */ 2106 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2107 2108 /* set "j" values into matrix */ 2109 maskcount = 1; 2110 for (j=0; j<nmask; j++) { 2111 a->j[jcount++] = masked[j]; 2112 mask[masked[j]] = maskcount++; 2113 } 2114 2115 /* set "a" values into matrix */ 2116 ishift = bs2*a->i[i]; 2117 for (j=0; j<bs; j++) { 2118 kmax = rowlengths[i*bs+j]; 2119 for (k=0; k<kmax; k++) { 2120 tmp = jj[nzcountb]/bs ; /* block col. index */ 2121 if (tmp >= i){ 2122 block = mask[tmp] - 1; 2123 point = jj[nzcountb] - bs*tmp; 2124 idx = ishift + bs2*block + j + bs*point; 2125 a->a[idx] = aa[nzcountb]; 2126 } 2127 nzcountb++; 2128 } 2129 } 2130 /* zero out the mask elements we set */ 2131 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2132 } 2133 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2134 2135 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2136 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2137 ierr = PetscFree(aa);CHKERRQ(ierr); 2138 ierr = PetscFree(jj);CHKERRQ(ierr); 2139 ierr = PetscFree(mask);CHKERRQ(ierr); 2140 2141 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2142 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2143 ierr = MatView_Private(B);CHKERRQ(ierr); 2144 *A = B; 2145 PetscFunctionReturn(0); 2146 } 2147 2148 2149 2150 #undef __FUNCT__ 2151 #define __FUNCT__ "MatRelax_SeqSBAIJ" 2152 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 2153 { 2154 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 2155 const MatScalar *aa=a->a,*v,*v1; 2156 PetscScalar *x,*b,*t,sum,d; 2157 MatScalar tmp; 2158 PetscErrorCode ierr; 2159 PetscInt m=a->mbs,bs=A->rmap->bs,j; 2160 const PetscInt *ai=a->i,*aj=a->j,*vj,*vj1; 2161 PetscInt nz,nz1,i; 2162 2163 PetscFunctionBegin; 2164 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 2165 2166 its = its*lits; 2167 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 2168 2169 if (bs > 1) SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 2170 2171 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 2172 if (xx != bb) { 2173 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2174 } else { 2175 b = x; 2176 } 2177 2178 if (!a->idiagvalid) { 2179 if (!a->idiag) { 2180 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 2181 } 2182 for (i=0; i<a->mbs; i++) a->idiag[i] = 1.0/a->a[a->i[i]]; 2183 a->idiagvalid = PETSC_TRUE; 2184 } 2185 2186 if (!a->relax_work) { 2187 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2188 } 2189 t = a->relax_work; 2190 2191 2192 if (flag & SOR_ZERO_INITIAL_GUESS) { 2193 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2194 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2195 2196 for (i=0; i<m; i++){ 2197 v = aa + ai[i] + 1; 2198 vj = aj + ai[i] + 1; 2199 nz = ai[i+1] - ai[i] - 1; 2200 x[i] = tmp = omega*t[i]*a->idiag[i]; 2201 for (j=0; j<nz; j++) { 2202 t[vj[j]] -= tmp*v[j]; 2203 } 2204 } 2205 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2206 } 2207 2208 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2209 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2210 t = b; 2211 } 2212 2213 for (i=m-1; i>=0; i--){ 2214 v = aa + ai[i] + 1; 2215 vj = aj + ai[i] + 1; 2216 nz = ai[i+1] - ai[i] - 1; 2217 sum = t[i]; 2218 PetscSparseDenseMinusDot(sum,x,v,vj,nz); 2219 x[i] = (1-omega)*x[i] + omega*sum*a->idiag[i]; 2220 } 2221 t = a->relax_work; 2222 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2223 } 2224 its--; 2225 } 2226 2227 while (its--) { 2228 /* 2229 forward sweep: 2230 for i=0,...,m-1: 2231 sum[i] = (b[i] - U(i,:)x )/d[i]; 2232 x[i] = (1-omega)x[i] + omega*sum[i]; 2233 b = b - x[i]*U^T(i,:); 2234 2235 */ 2236 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2237 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2238 2239 for (i=0; i<m; i++){ 2240 d = *(aa + ai[i]); /* diag[i] */ 2241 v = aa + ai[i] + 1; v1=v; 2242 vj = aj + ai[i] + 1; vj1=vj; 2243 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2244 sum = t[i]; 2245 ierr = PetscLogFlops(4.0*nz-2);CHKERRQ(ierr); 2246 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2247 x[i] = (1-omega)*x[i] + omega*sum/d; 2248 while (nz--) t[*vj++] -= x[i]*(*v++); 2249 } 2250 } 2251 2252 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2253 /* 2254 backward sweep: 2255 b = b - x[i]*U^T(i,:), i=0,...,n-2 2256 for i=m-1,...,0: 2257 sum[i] = (b[i] - U(i,:)x )/d[i]; 2258 x[i] = (1-omega)x[i] + omega*sum[i]; 2259 */ 2260 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2261 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2262 2263 for (i=0; i<m-1; i++){ /* update rhs */ 2264 v = aa + ai[i] + 1; 2265 vj = aj + ai[i] + 1; 2266 nz = ai[i+1] - ai[i] - 1; 2267 ierr = PetscLogFlops(2.0*nz-1);CHKERRQ(ierr); 2268 while (nz--) t[*vj++] -= x[i]*(*v++); 2269 } 2270 for (i=m-1; i>=0; i--){ 2271 d = *(aa + ai[i]); 2272 v = aa + ai[i] + 1; 2273 vj = aj + ai[i] + 1; 2274 nz = ai[i+1] - ai[i] - 1; 2275 ierr = PetscLogFlops(2.0*nz-1);CHKERRQ(ierr); 2276 sum = t[i]; 2277 while (nz--) sum -= x[*vj++]*(*v++); 2278 x[i] = (1-omega)*x[i] + omega*sum/d; 2279 } 2280 } 2281 } 2282 2283 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2284 if (bb != xx) { 2285 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2286 } 2287 PetscFunctionReturn(0); 2288 } 2289 2290 #undef __FUNCT__ 2291 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2292 /*@ 2293 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2294 (upper triangular entries in CSR format) provided by the user. 2295 2296 Collective on MPI_Comm 2297 2298 Input Parameters: 2299 + comm - must be an MPI communicator of size 1 2300 . bs - size of block 2301 . m - number of rows 2302 . n - number of columns 2303 . i - row indices 2304 . j - column indices 2305 - a - matrix values 2306 2307 Output Parameter: 2308 . mat - the matrix 2309 2310 Level: intermediate 2311 2312 Notes: 2313 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2314 once the matrix is destroyed 2315 2316 You cannot set new nonzero locations into this matrix, that will generate an error. 2317 2318 The i and j indices are 0 based 2319 2320 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2321 2322 @*/ 2323 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2324 { 2325 PetscErrorCode ierr; 2326 PetscInt ii; 2327 Mat_SeqSBAIJ *sbaij; 2328 2329 PetscFunctionBegin; 2330 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2331 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2332 2333 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2334 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2335 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2336 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2337 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2338 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2339 2340 sbaij->i = i; 2341 sbaij->j = j; 2342 sbaij->a = a; 2343 sbaij->singlemalloc = PETSC_FALSE; 2344 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2345 sbaij->free_a = PETSC_FALSE; 2346 sbaij->free_ij = PETSC_FALSE; 2347 2348 for (ii=0; ii<m; ii++) { 2349 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2350 #if defined(PETSC_USE_DEBUG) 2351 if (i[ii+1] - i[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %d length = %d",ii,i[ii+1] - i[ii]); 2352 #endif 2353 } 2354 #if defined(PETSC_USE_DEBUG) 2355 for (ii=0; ii<sbaij->i[m]; ii++) { 2356 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2357 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2358 } 2359 #endif 2360 2361 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2362 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2363 PetscFunctionReturn(0); 2364 } 2365 2366 2367 2368 2369 2370