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