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 /* 703 This is not yet used 704 */ 705 #undef __FUNCT__ 706 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_Inode" 707 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_Inode(Mat A) 708 { 709 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 710 PetscErrorCode ierr; 711 const PetscInt *ai = a->i, *aj = a->j,*cols; 712 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 713 PetscTruth flag; 714 715 PetscFunctionBegin; 716 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 717 while (i < m){ 718 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 719 /* Limits the number of elements in a node to 'a->inode.limit' */ 720 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 721 nzy = ai[j+1] - ai[j]; 722 if (nzy != (nzx - j + i)) break; 723 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 724 if (!flag) break; 725 } 726 ns[node_count++] = blk_size; 727 i = j; 728 } 729 if (!a->inode.size && m && node_count > .9*m) { 730 ierr = PetscFree(ns);CHKERRQ(ierr); 731 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 732 } else { 733 a->inode.node_count = node_count; 734 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 735 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt)); 736 ierr = PetscFree(ns);CHKERRQ(ierr); 737 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 738 739 /* count collections of adjacent columns in each inode */ 740 row = 0; 741 cnt = 0; 742 for (i=0; i<node_count; i++) { 743 cols = aj + ai[row] + a->inode.size[i]; 744 nz = ai[row+1] - ai[row] - a->inode.size[i]; 745 for (j=1; j<nz; j++) { 746 if (cols[j] != cols[j-1]+1) { 747 cnt++; 748 } 749 } 750 cnt++; 751 row += a->inode.size[i]; 752 } 753 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 754 cnt = 0; 755 row = 0; 756 for (i=0; i<node_count; i++) { 757 cols = aj + ai[row] + a->inode.size[i]; 758 CHKMEMQ; 759 counts[2*cnt] = cols[0]; 760 CHKMEMQ; 761 nz = ai[row+1] - ai[row] - a->inode.size[i]; 762 cnt2 = 1; 763 for (j=1; j<nz; j++) { 764 if (cols[j] != cols[j-1]+1) { 765 CHKMEMQ; 766 counts[2*(cnt++)+1] = cnt2; 767 counts[2*cnt] = cols[j]; 768 CHKMEMQ; 769 cnt2 = 1; 770 } else cnt2++; 771 } 772 CHKMEMQ; 773 counts[2*(cnt++)+1] = cnt2; 774 CHKMEMQ; 775 row += a->inode.size[i]; 776 } 777 ierr = PetscIntView(2*cnt,counts,0); 778 } 779 780 781 PetscFunctionReturn(0); 782 } 783 784 #undef __FUNCT__ 785 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 786 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 787 { 788 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 789 PetscErrorCode ierr; 790 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 791 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 792 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 793 MatScalar *aa = a->a,*ap; 794 795 PetscFunctionBegin; 796 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 797 798 if (m) rmax = ailen[0]; 799 for (i=1; i<mbs; i++) { 800 /* move each row back by the amount of empty slots (fshift) before it*/ 801 fshift += imax[i-1] - ailen[i-1]; 802 rmax = PetscMax(rmax,ailen[i]); 803 if (fshift) { 804 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 805 N = ailen[i]; 806 for (j=0; j<N; j++) { 807 ip[j-fshift] = ip[j]; 808 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 809 } 810 } 811 ai[i] = ai[i-1] + ailen[i-1]; 812 } 813 if (mbs) { 814 fshift += imax[mbs-1] - ailen[mbs-1]; 815 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 816 } 817 /* reset ilen and imax for each row */ 818 for (i=0; i<mbs; i++) { 819 ailen[i] = imax[i] = ai[i+1] - ai[i]; 820 } 821 a->nz = ai[mbs]; 822 823 /* diagonals may have moved, reset it */ 824 if (a->diag) { 825 ierr = PetscMemcpy(a->diag,ai,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 826 } 827 if (fshift && a->nounused == -1) { 828 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); 829 } 830 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); 831 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 832 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 833 a->reallocs = 0; 834 A->info.nz_unneeded = (PetscReal)fshift*bs2; 835 a->idiagvalid = PETSC_FALSE; 836 /* if (bs2 == 1) { 837 ierr = MatAssemblyEnd_SeqSBAIJ_Inode(A);CHKERRQ(ierr); 838 } */ 839 PetscFunctionReturn(0); 840 } 841 842 /* 843 This function returns an array of flags which indicate the locations of contiguous 844 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 845 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 846 Assume: sizes should be long enough to hold all the values. 847 */ 848 #undef __FUNCT__ 849 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 850 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 851 { 852 PetscInt i,j,k,row; 853 PetscTruth flg; 854 855 PetscFunctionBegin; 856 for (i=0,j=0; i<n; j++) { 857 row = idx[i]; 858 if (row%bs!=0) { /* Not the begining of a block */ 859 sizes[j] = 1; 860 i++; 861 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 862 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 863 i++; 864 } else { /* Begining of the block, so check if the complete block exists */ 865 flg = PETSC_TRUE; 866 for (k=1; k<bs; k++) { 867 if (row+k != idx[i+k]) { /* break in the block */ 868 flg = PETSC_FALSE; 869 break; 870 } 871 } 872 if (flg) { /* No break in the bs */ 873 sizes[j] = bs; 874 i+= bs; 875 } else { 876 sizes[j] = 1; 877 i++; 878 } 879 } 880 } 881 *bs_max = j; 882 PetscFunctionReturn(0); 883 } 884 885 886 /* Only add/insert a(i,j) with i<=j (blocks). 887 Any a(i,j) with i>j input by user is ingored. 888 */ 889 890 #undef __FUNCT__ 891 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 892 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 893 { 894 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 895 PetscErrorCode ierr; 896 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 897 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 898 PetscInt *aj=a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 899 PetscInt ridx,cidx,bs2=a->bs2; 900 MatScalar *ap,value,*aa=a->a,*bap; 901 902 PetscFunctionBegin; 903 for (k=0; k<m; k++) { /* loop over added rows */ 904 row = im[k]; /* row number */ 905 brow = row/bs; /* block row number */ 906 if (row < 0) continue; 907 #if defined(PETSC_USE_DEBUG) 908 if (row >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,A->rmap->N-1); 909 #endif 910 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 911 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 912 rmax = imax[brow]; /* maximum space allocated for this row */ 913 nrow = ailen[brow]; /* actual length of this row */ 914 low = 0; 915 916 for (l=0; l<n; l++) { /* loop over added columns */ 917 if (in[l] < 0) continue; 918 #if defined(PETSC_USE_DEBUG) 919 if (in[l] >= A->rmap->N) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",in[l],A->rmap->N-1); 920 #endif 921 col = in[l]; 922 bcol = col/bs; /* block col number */ 923 924 if (brow > bcol) { 925 if (a->ignore_ltriangular){ 926 continue; /* ignore lower triangular values */ 927 } else { 928 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)"); 929 } 930 } 931 932 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 933 if ((brow==bcol && ridx<=cidx) || (brow<bcol)){ 934 /* element value a(k,l) */ 935 if (roworiented) { 936 value = v[l + k*n]; 937 } else { 938 value = v[k + l*m]; 939 } 940 941 /* move pointer bap to a(k,l) quickly and add/insert value */ 942 if (col <= lastcol) low = 0; high = nrow; 943 lastcol = col; 944 while (high-low > 7) { 945 t = (low+high)/2; 946 if (rp[t] > bcol) high = t; 947 else low = t; 948 } 949 for (i=low; i<high; i++) { 950 if (rp[i] > bcol) break; 951 if (rp[i] == bcol) { 952 bap = ap + bs2*i + bs*cidx + ridx; 953 if (is == ADD_VALUES) *bap += value; 954 else *bap = value; 955 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 956 if (brow == bcol && ridx < cidx){ 957 bap = ap + bs2*i + bs*ridx + cidx; 958 if (is == ADD_VALUES) *bap += value; 959 else *bap = value; 960 } 961 goto noinsert1; 962 } 963 } 964 965 if (nonew == 1) goto noinsert1; 966 if (nonew == -1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 967 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 968 969 N = nrow++ - 1; high++; 970 /* shift up all the later entries in this row */ 971 for (ii=N; ii>=i; ii--) { 972 rp[ii+1] = rp[ii]; 973 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 974 } 975 if (N>=i) { 976 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 977 } 978 rp[i] = bcol; 979 ap[bs2*i + bs*cidx + ridx] = value; 980 noinsert1:; 981 low = i; 982 } 983 } /* end of loop over added columns */ 984 ailen[brow] = nrow; 985 } /* end of loop over added rows */ 986 PetscFunctionReturn(0); 987 } 988 989 #undef __FUNCT__ 990 #define __FUNCT__ "MatICCFactor_SeqSBAIJ" 991 PetscErrorCode MatICCFactor_SeqSBAIJ(Mat inA,IS row,const MatFactorInfo *info) 992 { 993 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)inA->data; 994 Mat outA; 995 PetscErrorCode ierr; 996 PetscTruth row_identity; 997 998 PetscFunctionBegin; 999 if (info->levels != 0) SETERRQ(PETSC_ERR_SUP,"Only levels=0 is supported for in-place icc"); 1000 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 1001 if (!row_identity) SETERRQ(PETSC_ERR_SUP,"Matrix reordering is not supported"); 1002 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()! */ 1003 1004 outA = inA; 1005 inA->factor = MAT_FACTOR_ICC; 1006 1007 ierr = MatMarkDiagonal_SeqSBAIJ(inA);CHKERRQ(ierr); 1008 ierr = MatSeqSBAIJSetNumericFactorization(inA,row_identity);CHKERRQ(ierr); 1009 1010 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1011 if (a->row) { ierr = ISDestroy(a->row);CHKERRQ(ierr); } 1012 a->row = row; 1013 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 1014 if (a->col) { ierr = ISDestroy(a->col);CHKERRQ(ierr); } 1015 a->col = row; 1016 1017 /* Create the invert permutation so that it can be used in MatCholeskyFactorNumeric() */ 1018 if (a->icol) {ierr = ISInvertPermutation(row,PETSC_DECIDE, &a->icol);CHKERRQ(ierr);} 1019 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 1020 1021 if (!a->solve_work) { 1022 ierr = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 1023 ierr = PetscLogObjectMemory(inA,(inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar));CHKERRQ(ierr); 1024 } 1025 1026 ierr = MatCholeskyFactorNumeric(outA,inA,info);CHKERRQ(ierr); 1027 PetscFunctionReturn(0); 1028 } 1029 1030 EXTERN_C_BEGIN 1031 #undef __FUNCT__ 1032 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices_SeqSBAIJ" 1033 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices_SeqSBAIJ(Mat mat,PetscInt *indices) 1034 { 1035 Mat_SeqSBAIJ *baij = (Mat_SeqSBAIJ *)mat->data; 1036 PetscInt i,nz,n; 1037 1038 PetscFunctionBegin; 1039 nz = baij->maxnz; 1040 n = mat->cmap->n; 1041 for (i=0; i<nz; i++) { 1042 baij->j[i] = indices[i]; 1043 } 1044 baij->nz = nz; 1045 for (i=0; i<n; i++) { 1046 baij->ilen[i] = baij->imax[i]; 1047 } 1048 PetscFunctionReturn(0); 1049 } 1050 EXTERN_C_END 1051 1052 #undef __FUNCT__ 1053 #define __FUNCT__ "MatSeqSBAIJSetColumnIndices" 1054 /*@ 1055 MatSeqSBAIJSetColumnIndices - Set the column indices for all the rows 1056 in the matrix. 1057 1058 Input Parameters: 1059 + mat - the SeqSBAIJ matrix 1060 - indices - the column indices 1061 1062 Level: advanced 1063 1064 Notes: 1065 This can be called if you have precomputed the nonzero structure of the 1066 matrix and want to provide it to the matrix object to improve the performance 1067 of the MatSetValues() operation. 1068 1069 You MUST have set the correct numbers of nonzeros per row in the call to 1070 MatCreateSeqSBAIJ(), and the columns indices MUST be sorted. 1071 1072 MUST be called before any calls to MatSetValues() 1073 1074 .seealso: MatCreateSeqSBAIJ 1075 @*/ 1076 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetColumnIndices(Mat mat,PetscInt *indices) 1077 { 1078 PetscErrorCode ierr,(*f)(Mat,PetscInt *); 1079 1080 PetscFunctionBegin; 1081 PetscValidHeaderSpecific(mat,MAT_COOKIE,1); 1082 PetscValidPointer(indices,2); 1083 ierr = PetscObjectQueryFunction((PetscObject)mat,"MatSeqSBAIJSetColumnIndices_C",(void (**)(void))&f);CHKERRQ(ierr); 1084 if (f) { 1085 ierr = (*f)(mat,indices);CHKERRQ(ierr); 1086 } else { 1087 SETERRQ(PETSC_ERR_SUP,"Wrong type of matrix to set column indices"); 1088 } 1089 PetscFunctionReturn(0); 1090 } 1091 1092 #undef __FUNCT__ 1093 #define __FUNCT__ "MatCopy_SeqSBAIJ" 1094 PetscErrorCode MatCopy_SeqSBAIJ(Mat A,Mat B,MatStructure str) 1095 { 1096 PetscErrorCode ierr; 1097 1098 PetscFunctionBegin; 1099 /* If the two matrices have the same copy implementation, use fast copy. */ 1100 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 1101 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1102 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1103 1104 if (a->i[A->rmap->N] != b->i[B->rmap->N]) { 1105 SETERRQ(PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different"); 1106 } 1107 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->N])*sizeof(PetscScalar));CHKERRQ(ierr); 1108 } else { 1109 ierr = MatGetRowUpperTriangular(A);CHKERRQ(ierr); 1110 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 1111 ierr = MatRestoreRowUpperTriangular(A);CHKERRQ(ierr); 1112 } 1113 PetscFunctionReturn(0); 1114 } 1115 1116 #undef __FUNCT__ 1117 #define __FUNCT__ "MatSetUpPreallocation_SeqSBAIJ" 1118 PetscErrorCode MatSetUpPreallocation_SeqSBAIJ(Mat A) 1119 { 1120 PetscErrorCode ierr; 1121 1122 PetscFunctionBegin; 1123 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(A,-PetscMax(A->rmap->bs,1),PETSC_DEFAULT,0);CHKERRQ(ierr); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 #undef __FUNCT__ 1128 #define __FUNCT__ "MatGetArray_SeqSBAIJ" 1129 PetscErrorCode MatGetArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1130 { 1131 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1132 PetscFunctionBegin; 1133 *array = a->a; 1134 PetscFunctionReturn(0); 1135 } 1136 1137 #undef __FUNCT__ 1138 #define __FUNCT__ "MatRestoreArray_SeqSBAIJ" 1139 PetscErrorCode MatRestoreArray_SeqSBAIJ(Mat A,PetscScalar *array[]) 1140 { 1141 PetscFunctionBegin; 1142 PetscFunctionReturn(0); 1143 } 1144 1145 #include "petscblaslapack.h" 1146 #undef __FUNCT__ 1147 #define __FUNCT__ "MatAXPY_SeqSBAIJ" 1148 PetscErrorCode MatAXPY_SeqSBAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 1149 { 1150 Mat_SeqSBAIJ *x=(Mat_SeqSBAIJ *)X->data, *y=(Mat_SeqSBAIJ *)Y->data; 1151 PetscErrorCode ierr; 1152 PetscInt i,bs=Y->rmap->bs,bs2,j; 1153 PetscBLASInt one = 1,bnz = PetscBLASIntCast(x->nz); 1154 1155 PetscFunctionBegin; 1156 if (str == SAME_NONZERO_PATTERN) { 1157 PetscScalar alpha = a; 1158 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 1159 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 1160 if (y->xtoy && y->XtoY != X) { 1161 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 1162 ierr = MatDestroy(y->XtoY);CHKERRQ(ierr); 1163 } 1164 if (!y->xtoy) { /* get xtoy */ 1165 ierr = MatAXPYGetxtoy_Private(x->mbs,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 1166 y->XtoY = X; 1167 } 1168 bs2 = bs*bs; 1169 for (i=0; i<x->nz; i++) { 1170 j = 0; 1171 while (j < bs2){ 1172 y->a[bs2*y->xtoy[i]+j] += a*(x->a[bs2*i+j]); 1173 j++; 1174 } 1175 } 1176 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); 1177 } else { 1178 ierr = MatGetRowUpperTriangular(X);CHKERRQ(ierr); 1179 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 1180 ierr = MatRestoreRowUpperTriangular(X);CHKERRQ(ierr); 1181 } 1182 PetscFunctionReturn(0); 1183 } 1184 1185 #undef __FUNCT__ 1186 #define __FUNCT__ "MatIsSymmetric_SeqSBAIJ" 1187 PetscErrorCode MatIsSymmetric_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1188 { 1189 PetscFunctionBegin; 1190 *flg = PETSC_TRUE; 1191 PetscFunctionReturn(0); 1192 } 1193 1194 #undef __FUNCT__ 1195 #define __FUNCT__ "MatIsStructurallySymmetric_SeqSBAIJ" 1196 PetscErrorCode MatIsStructurallySymmetric_SeqSBAIJ(Mat A,PetscTruth *flg) 1197 { 1198 PetscFunctionBegin; 1199 *flg = PETSC_TRUE; 1200 PetscFunctionReturn(0); 1201 } 1202 1203 #undef __FUNCT__ 1204 #define __FUNCT__ "MatIsHermitian_SeqSBAIJ" 1205 PetscErrorCode MatIsHermitian_SeqSBAIJ(Mat A,PetscReal tol,PetscTruth *flg) 1206 { 1207 PetscFunctionBegin; 1208 *flg = PETSC_FALSE; 1209 PetscFunctionReturn(0); 1210 } 1211 1212 #undef __FUNCT__ 1213 #define __FUNCT__ "MatRealPart_SeqSBAIJ" 1214 PetscErrorCode MatRealPart_SeqSBAIJ(Mat A) 1215 { 1216 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1217 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1218 MatScalar *aa = a->a; 1219 1220 PetscFunctionBegin; 1221 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1222 PetscFunctionReturn(0); 1223 } 1224 1225 #undef __FUNCT__ 1226 #define __FUNCT__ "MatImaginaryPart_SeqSBAIJ" 1227 PetscErrorCode MatImaginaryPart_SeqSBAIJ(Mat A) 1228 { 1229 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 1230 PetscInt i,nz = a->bs2*a->i[a->mbs]; 1231 MatScalar *aa = a->a; 1232 1233 PetscFunctionBegin; 1234 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1235 PetscFunctionReturn(0); 1236 } 1237 1238 /* -------------------------------------------------------------------*/ 1239 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1240 MatGetRow_SeqSBAIJ, 1241 MatRestoreRow_SeqSBAIJ, 1242 MatMult_SeqSBAIJ_N, 1243 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1244 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1245 MatMultAdd_SeqSBAIJ_N, 1246 0, 1247 0, 1248 0, 1249 /*10*/ 0, 1250 0, 1251 MatCholeskyFactor_SeqSBAIJ, 1252 MatRelax_SeqSBAIJ, 1253 MatTranspose_SeqSBAIJ, 1254 /*15*/ MatGetInfo_SeqSBAIJ, 1255 MatEqual_SeqSBAIJ, 1256 MatGetDiagonal_SeqSBAIJ, 1257 MatDiagonalScale_SeqSBAIJ, 1258 MatNorm_SeqSBAIJ, 1259 /*20*/ 0, 1260 MatAssemblyEnd_SeqSBAIJ, 1261 MatSetOption_SeqSBAIJ, 1262 MatZeroEntries_SeqSBAIJ, 1263 /*24*/ 0, 1264 0, 1265 0, 1266 0, 1267 0, 1268 /*29*/ MatSetUpPreallocation_SeqSBAIJ, 1269 0, 1270 0, 1271 MatGetArray_SeqSBAIJ, 1272 MatRestoreArray_SeqSBAIJ, 1273 /*34*/ MatDuplicate_SeqSBAIJ, 1274 0, 1275 0, 1276 0, 1277 MatICCFactor_SeqSBAIJ, 1278 /*39*/ MatAXPY_SeqSBAIJ, 1279 MatGetSubMatrices_SeqSBAIJ, 1280 MatIncreaseOverlap_SeqSBAIJ, 1281 MatGetValues_SeqSBAIJ, 1282 MatCopy_SeqSBAIJ, 1283 /*44*/ 0, 1284 MatScale_SeqSBAIJ, 1285 0, 1286 0, 1287 0, 1288 /*49*/ 0, 1289 MatGetRowIJ_SeqSBAIJ, 1290 MatRestoreRowIJ_SeqSBAIJ, 1291 0, 1292 0, 1293 /*54*/ 0, 1294 0, 1295 0, 1296 0, 1297 MatSetValuesBlocked_SeqSBAIJ, 1298 /*59*/ MatGetSubMatrix_SeqSBAIJ, 1299 0, 1300 0, 1301 0, 1302 0, 1303 /*64*/ 0, 1304 0, 1305 0, 1306 0, 1307 0, 1308 /*69*/ MatGetRowMaxAbs_SeqSBAIJ, 1309 0, 1310 0, 1311 0, 1312 0, 1313 /*74*/ 0, 1314 0, 1315 0, 1316 0, 1317 0, 1318 /*79*/ 0, 1319 0, 1320 0, 1321 #if !defined(PETSC_USE_COMPLEX) 1322 MatGetInertia_SeqSBAIJ, 1323 #else 1324 0, 1325 #endif 1326 MatLoad_SeqSBAIJ, 1327 /*84*/ MatIsSymmetric_SeqSBAIJ, 1328 MatIsHermitian_SeqSBAIJ, 1329 MatIsStructurallySymmetric_SeqSBAIJ, 1330 0, 1331 0, 1332 /*89*/ 0, 1333 0, 1334 0, 1335 0, 1336 0, 1337 /*94*/ 0, 1338 0, 1339 0, 1340 0, 1341 0, 1342 /*99*/ 0, 1343 0, 1344 0, 1345 0, 1346 0, 1347 /*104*/0, 1348 MatRealPart_SeqSBAIJ, 1349 MatImaginaryPart_SeqSBAIJ, 1350 MatGetRowUpperTriangular_SeqSBAIJ, 1351 MatRestoreRowUpperTriangular_SeqSBAIJ, 1352 /*109*/0, 1353 0, 1354 0, 1355 0, 1356 MatMissingDiagonal_SeqSBAIJ 1357 }; 1358 1359 EXTERN_C_BEGIN 1360 #undef __FUNCT__ 1361 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1362 PetscErrorCode PETSCMAT_DLLEXPORT MatStoreValues_SeqSBAIJ(Mat mat) 1363 { 1364 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1365 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1366 PetscErrorCode ierr; 1367 1368 PetscFunctionBegin; 1369 if (aij->nonew != 1) { 1370 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1371 } 1372 1373 /* allocate space for values if not already there */ 1374 if (!aij->saved_values) { 1375 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1376 } 1377 1378 /* copy values over */ 1379 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1380 PetscFunctionReturn(0); 1381 } 1382 EXTERN_C_END 1383 1384 EXTERN_C_BEGIN 1385 #undef __FUNCT__ 1386 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1387 PetscErrorCode PETSCMAT_DLLEXPORT MatRetrieveValues_SeqSBAIJ(Mat mat) 1388 { 1389 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ *)mat->data; 1390 PetscErrorCode ierr; 1391 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1392 1393 PetscFunctionBegin; 1394 if (aij->nonew != 1) { 1395 SETERRQ(PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1396 } 1397 if (!aij->saved_values) { 1398 SETERRQ(PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1399 } 1400 1401 /* copy values over */ 1402 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1403 PetscFunctionReturn(0); 1404 } 1405 EXTERN_C_END 1406 1407 EXTERN_C_BEGIN 1408 #undef __FUNCT__ 1409 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1410 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1411 { 1412 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1413 PetscErrorCode ierr; 1414 PetscInt i,mbs,bs2, newbs = PetscAbs(bs); 1415 PetscTruth skipallocation = PETSC_FALSE,flg = PETSC_FALSE; 1416 1417 PetscFunctionBegin; 1418 B->preallocated = PETSC_TRUE; 1419 if (bs < 0) { 1420 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for MPISBAIJ matrix","Mat");CHKERRQ(ierr); 1421 ierr = PetscOptionsInt("-mat_block_size","Set the blocksize used to store the matrix","MatSeqSBAIJSetPreallocation",newbs,&newbs,PETSC_NULL);CHKERRQ(ierr); 1422 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1423 bs = PetscAbs(bs); 1424 } 1425 if (nnz && newbs != bs) { 1426 SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot change blocksize from command line if setting nnz"); 1427 } 1428 bs = newbs; 1429 1430 ierr = PetscMapSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1431 ierr = PetscMapSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1432 ierr = PetscMapSetUp(B->rmap);CHKERRQ(ierr); 1433 ierr = PetscMapSetUp(B->cmap);CHKERRQ(ierr); 1434 1435 mbs = B->rmap->N/bs; 1436 bs2 = bs*bs; 1437 1438 if (mbs*bs != B->rmap->N) { 1439 SETERRQ(PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1440 } 1441 1442 if (nz == MAT_SKIP_ALLOCATION) { 1443 skipallocation = PETSC_TRUE; 1444 nz = 0; 1445 } 1446 1447 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1448 if (nz < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1449 if (nnz) { 1450 for (i=0; i<mbs; i++) { 1451 if (nnz[i] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %D value %D",i,nnz[i]); 1452 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); 1453 } 1454 } 1455 1456 B->ops->mult = MatMult_SeqSBAIJ_N; 1457 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1458 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1459 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1460 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1461 if (!flg) { 1462 switch (bs) { 1463 case 1: 1464 B->ops->mult = MatMult_SeqSBAIJ_1; 1465 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1466 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1467 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1468 break; 1469 case 2: 1470 B->ops->mult = MatMult_SeqSBAIJ_2; 1471 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1472 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1473 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1474 break; 1475 case 3: 1476 B->ops->mult = MatMult_SeqSBAIJ_3; 1477 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1478 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1479 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1480 break; 1481 case 4: 1482 B->ops->mult = MatMult_SeqSBAIJ_4; 1483 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1484 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1485 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1486 break; 1487 case 5: 1488 B->ops->mult = MatMult_SeqSBAIJ_5; 1489 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1490 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1491 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1492 break; 1493 case 6: 1494 B->ops->mult = MatMult_SeqSBAIJ_6; 1495 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1496 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1497 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1498 break; 1499 case 7: 1500 B->ops->mult = MatMult_SeqSBAIJ_7; 1501 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1502 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1503 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1504 break; 1505 } 1506 } 1507 1508 b->mbs = mbs; 1509 b->nbs = mbs; 1510 if (!skipallocation) { 1511 if (!b->imax) { 1512 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1513 ierr = PetscLogObjectMemory(B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1514 } 1515 if (!nnz) { 1516 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1517 else if (nz <= 0) nz = 1; 1518 for (i=0; i<mbs; i++) { 1519 b->imax[i] = nz; 1520 } 1521 nz = nz*mbs; /* total nz */ 1522 } else { 1523 nz = 0; 1524 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1525 } 1526 /* b->ilen will count nonzeros in each block row so far. */ 1527 for (i=0; i<mbs; i++) { b->ilen[i] = 0;} 1528 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1529 1530 /* allocate the matrix space */ 1531 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1532 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1533 ierr = PetscLogObjectMemory(B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1534 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1535 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1536 b->singlemalloc = PETSC_TRUE; 1537 1538 /* pointer to beginning of each row */ 1539 b->i[0] = 0; 1540 for (i=1; i<mbs+1; i++) { 1541 b->i[i] = b->i[i-1] + b->imax[i-1]; 1542 } 1543 b->free_a = PETSC_TRUE; 1544 b->free_ij = PETSC_TRUE; 1545 } else { 1546 b->free_a = PETSC_FALSE; 1547 b->free_ij = PETSC_FALSE; 1548 } 1549 1550 B->rmap->bs = bs; 1551 b->bs2 = bs2; 1552 b->nz = 0; 1553 b->maxnz = nz*bs2; 1554 1555 b->inew = 0; 1556 b->jnew = 0; 1557 b->anew = 0; 1558 b->a2anew = 0; 1559 b->permute = PETSC_FALSE; 1560 PetscFunctionReturn(0); 1561 } 1562 EXTERN_C_END 1563 1564 #undef __FUNCT__ 1565 #define __FUNCT__ "MatSeqBAIJSetNumericFactorization" 1566 /* 1567 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1568 */ 1569 PetscErrorCode MatSeqSBAIJSetNumericFactorization(Mat B,PetscTruth natural) 1570 { 1571 PetscErrorCode ierr; 1572 PetscTruth flg = PETSC_FALSE; 1573 PetscInt bs = B->rmap->bs; 1574 1575 PetscFunctionBegin; 1576 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,PETSC_NULL);CHKERRQ(ierr); 1577 if (flg) bs = 8; 1578 1579 if (!natural) { 1580 switch (bs) { 1581 case 1: 1582 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1; 1583 break; 1584 case 2: 1585 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1586 break; 1587 case 3: 1588 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1589 break; 1590 case 4: 1591 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1592 break; 1593 case 5: 1594 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1595 break; 1596 case 6: 1597 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1598 break; 1599 case 7: 1600 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1601 break; 1602 default: 1603 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1604 break; 1605 } 1606 } else { 1607 switch (bs) { 1608 case 1: 1609 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering; 1610 break; 1611 case 2: 1612 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1613 break; 1614 case 3: 1615 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1616 break; 1617 case 4: 1618 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1619 break; 1620 case 5: 1621 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1622 break; 1623 case 6: 1624 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1625 break; 1626 case 7: 1627 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1628 break; 1629 default: 1630 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1631 break; 1632 } 1633 } 1634 PetscFunctionReturn(0); 1635 } 1636 1637 EXTERN_C_BEGIN 1638 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1639 EXTERN PetscErrorCode PETSCMAT_DLLEXPORT MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1640 EXTERN_C_END 1641 1642 1643 EXTERN_C_BEGIN 1644 #undef __FUNCT__ 1645 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1646 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1647 { 1648 PetscInt n = A->rmap->n; 1649 PetscErrorCode ierr; 1650 1651 PetscFunctionBegin; 1652 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1653 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1654 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1655 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1656 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1657 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1658 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1659 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1660 (*B)->factor = ftype; 1661 PetscFunctionReturn(0); 1662 } 1663 EXTERN_C_END 1664 1665 EXTERN_C_BEGIN 1666 #undef __FUNCT__ 1667 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1668 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1669 { 1670 PetscFunctionBegin; 1671 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1672 *flg = PETSC_TRUE; 1673 } else { 1674 *flg = PETSC_FALSE; 1675 } 1676 PetscFunctionReturn(0); 1677 } 1678 EXTERN_C_END 1679 1680 EXTERN_C_BEGIN 1681 #if defined(PETSC_HAVE_MUMPS) 1682 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1683 #endif 1684 #if defined(PETSC_HAVE_SPOOLES) 1685 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1686 #endif 1687 #if defined(PETSC_HAVE_PASTIX) 1688 extern PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1689 #endif 1690 EXTERN_C_END 1691 1692 /*MC 1693 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1694 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1695 1696 Options Database Keys: 1697 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1698 1699 Level: beginner 1700 1701 .seealso: MatCreateSeqSBAIJ 1702 M*/ 1703 1704 EXTERN_C_BEGIN 1705 #undef __FUNCT__ 1706 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1707 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1708 { 1709 Mat_SeqSBAIJ *b; 1710 PetscErrorCode ierr; 1711 PetscMPIInt size; 1712 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1713 1714 PetscFunctionBegin; 1715 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1716 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1717 1718 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1719 B->data = (void*)b; 1720 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1721 B->ops->destroy = MatDestroy_SeqSBAIJ; 1722 B->ops->view = MatView_SeqSBAIJ; 1723 B->mapping = 0; 1724 b->row = 0; 1725 b->icol = 0; 1726 b->reallocs = 0; 1727 b->saved_values = 0; 1728 b->inode.limit = 5; 1729 b->inode.max_limit = 5; 1730 1731 b->roworiented = PETSC_TRUE; 1732 b->nonew = 0; 1733 b->diag = 0; 1734 b->solve_work = 0; 1735 b->mult_work = 0; 1736 B->spptr = 0; 1737 b->keepnonzeropattern = PETSC_FALSE; 1738 b->xtoy = 0; 1739 b->XtoY = 0; 1740 1741 b->inew = 0; 1742 b->jnew = 0; 1743 b->anew = 0; 1744 b->a2anew = 0; 1745 b->permute = PETSC_FALSE; 1746 1747 b->ignore_ltriangular = PETSC_FALSE; 1748 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1749 1750 b->getrow_utriangular = PETSC_FALSE; 1751 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1752 1753 #if defined(PETSC_HAVE_PASTIX) 1754 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_pastix_C", 1755 "MatGetFactor_seqsbaij_pastix", 1756 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1757 #endif 1758 #if defined(PETSC_HAVE_SPOOLES) 1759 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_spooles_C", 1760 "MatGetFactor_seqsbaij_spooles", 1761 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1762 #endif 1763 #if defined(PETSC_HAVE_MUMPS) 1764 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_mumps_C", 1765 "MatGetFactor_seqsbaij_mumps", 1766 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1767 #endif 1768 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_seqsbaij_petsc_C", 1769 "MatGetFactorAvailable_seqsbaij_petsc", 1770 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1771 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_petsc_C", 1772 "MatGetFactor_seqsbaij_petsc", 1773 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1774 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1775 "MatStoreValues_SeqSBAIJ", 1776 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1777 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1778 "MatRetrieveValues_SeqSBAIJ", 1779 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1780 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1781 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1782 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1783 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1784 "MatConvert_SeqSBAIJ_SeqAIJ", 1785 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1786 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1787 "MatConvert_SeqSBAIJ_SeqBAIJ", 1788 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1789 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1790 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1791 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1792 1793 B->symmetric = PETSC_TRUE; 1794 B->structurally_symmetric = PETSC_TRUE; 1795 B->symmetric_set = PETSC_TRUE; 1796 B->structurally_symmetric_set = PETSC_TRUE; 1797 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1798 1799 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1800 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1801 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1802 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1803 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1804 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); 1805 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1806 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1807 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1808 1809 PetscFunctionReturn(0); 1810 } 1811 EXTERN_C_END 1812 1813 #undef __FUNCT__ 1814 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1815 /*@C 1816 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1817 compressed row) format. For good matrix assembly performance the 1818 user should preallocate the matrix storage by setting the parameter nz 1819 (or the array nnz). By setting these parameters accurately, performance 1820 during matrix assembly can be increased by more than a factor of 50. 1821 1822 Collective on Mat 1823 1824 Input Parameters: 1825 + A - the symmetric matrix 1826 . bs - size of block 1827 . nz - number of block nonzeros per block row (same for all rows) 1828 - nnz - array containing the number of block nonzeros in the upper triangular plus 1829 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1830 1831 Options Database Keys: 1832 . -mat_no_unroll - uses code that does not unroll the loops in the 1833 block calculations (much slower) 1834 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1835 1836 Level: intermediate 1837 1838 Notes: 1839 Specify the preallocated storage with either nz or nnz (not both). 1840 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1841 allocation. For additional details, see the users manual chapter on 1842 matrices. 1843 1844 You can call MatGetInfo() to get information on how effective the preallocation was; 1845 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1846 You can also run with the option -info and look for messages with the string 1847 malloc in them to see if additional memory allocation was needed. 1848 1849 If the nnz parameter is given then the nz parameter is ignored 1850 1851 1852 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1853 @*/ 1854 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1855 { 1856 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1857 1858 PetscFunctionBegin; 1859 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1860 if (f) { 1861 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1862 } 1863 PetscFunctionReturn(0); 1864 } 1865 1866 #undef __FUNCT__ 1867 #define __FUNCT__ "MatCreateSeqSBAIJ" 1868 /*@C 1869 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1870 compressed row) format. For good matrix assembly performance the 1871 user should preallocate the matrix storage by setting the parameter nz 1872 (or the array nnz). By setting these parameters accurately, performance 1873 during matrix assembly can be increased by more than a factor of 50. 1874 1875 Collective on MPI_Comm 1876 1877 Input Parameters: 1878 + comm - MPI communicator, set to PETSC_COMM_SELF 1879 . bs - size of block 1880 . m - number of rows, or number of columns 1881 . nz - number of block nonzeros per block row (same for all rows) 1882 - nnz - array containing the number of block nonzeros in the upper triangular plus 1883 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1884 1885 Output Parameter: 1886 . A - the symmetric matrix 1887 1888 Options Database Keys: 1889 . -mat_no_unroll - uses code that does not unroll the loops in the 1890 block calculations (much slower) 1891 . -mat_block_size - size of the blocks to use 1892 1893 Level: intermediate 1894 1895 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1896 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1897 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1898 1899 Notes: 1900 The number of rows and columns must be divisible by blocksize. 1901 This matrix type does not support complex Hermitian operation. 1902 1903 Specify the preallocated storage with either nz or nnz (not both). 1904 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1905 allocation. For additional details, see the users manual chapter on 1906 matrices. 1907 1908 If the nnz parameter is given then the nz parameter is ignored 1909 1910 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1911 @*/ 1912 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1913 { 1914 PetscErrorCode ierr; 1915 1916 PetscFunctionBegin; 1917 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1918 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1919 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1920 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1921 PetscFunctionReturn(0); 1922 } 1923 1924 #undef __FUNCT__ 1925 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1926 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1927 { 1928 Mat C; 1929 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1930 PetscErrorCode ierr; 1931 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1932 1933 PetscFunctionBegin; 1934 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1935 1936 *B = 0; 1937 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1938 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1939 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1940 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1941 c = (Mat_SeqSBAIJ*)C->data; 1942 1943 C->preallocated = PETSC_TRUE; 1944 C->factor = A->factor; 1945 c->row = 0; 1946 c->icol = 0; 1947 c->saved_values = 0; 1948 c->keepnonzeropattern = a->keepnonzeropattern; 1949 C->assembled = PETSC_TRUE; 1950 1951 ierr = PetscMapCopy(((PetscObject)A)->comm,A->rmap,C->rmap);CHKERRQ(ierr); 1952 ierr = PetscMapCopy(((PetscObject)A)->comm,A->cmap,C->cmap);CHKERRQ(ierr); 1953 c->bs2 = a->bs2; 1954 c->mbs = a->mbs; 1955 c->nbs = a->nbs; 1956 1957 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1958 for (i=0; i<mbs; i++) { 1959 c->imax[i] = a->imax[i]; 1960 c->ilen[i] = a->ilen[i]; 1961 } 1962 1963 /* allocate the matrix space */ 1964 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1965 c->singlemalloc = PETSC_TRUE; 1966 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1967 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1968 if (mbs > 0) { 1969 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1970 if (cpvalues == MAT_COPY_VALUES) { 1971 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1972 } else { 1973 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1974 } 1975 } 1976 1977 c->roworiented = a->roworiented; 1978 c->nonew = a->nonew; 1979 1980 if (a->diag) { 1981 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1982 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1983 for (i=0; i<mbs; i++) { 1984 c->diag[i] = a->diag[i]; 1985 } 1986 } else c->diag = 0; 1987 c->nz = a->nz; 1988 c->maxnz = a->maxnz; 1989 c->solve_work = 0; 1990 c->mult_work = 0; 1991 c->free_a = PETSC_TRUE; 1992 c->free_ij = PETSC_TRUE; 1993 *B = C; 1994 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 1995 PetscFunctionReturn(0); 1996 } 1997 1998 #undef __FUNCT__ 1999 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2000 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 2001 { 2002 Mat_SeqSBAIJ *a; 2003 Mat B; 2004 PetscErrorCode ierr; 2005 int fd; 2006 PetscMPIInt size; 2007 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2008 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2009 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2010 PetscInt *masked,nmask,tmp,bs2,ishift; 2011 PetscScalar *aa; 2012 MPI_Comm comm = ((PetscObject)viewer)->comm; 2013 2014 PetscFunctionBegin; 2015 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2016 bs2 = bs*bs; 2017 2018 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2019 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 2020 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2021 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2022 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2023 M = header[1]; N = header[2]; nz = header[3]; 2024 2025 if (header[3] < 0) { 2026 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2027 } 2028 2029 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 2030 2031 /* 2032 This code adds extra rows to make sure the number of rows is 2033 divisible by the blocksize 2034 */ 2035 mbs = M/bs; 2036 extra_rows = bs - M + bs*(mbs); 2037 if (extra_rows == bs) extra_rows = 0; 2038 else mbs++; 2039 if (extra_rows) { 2040 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2041 } 2042 2043 /* read in row lengths */ 2044 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2045 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2046 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2047 2048 /* read in column indices */ 2049 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2050 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2051 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2052 2053 /* loop over row lengths determining block row lengths */ 2054 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2055 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2056 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 2057 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2058 masked = mask + mbs; 2059 rowcount = 0; nzcount = 0; 2060 for (i=0; i<mbs; i++) { 2061 nmask = 0; 2062 for (j=0; j<bs; j++) { 2063 kmax = rowlengths[rowcount]; 2064 for (k=0; k<kmax; k++) { 2065 tmp = jj[nzcount++]/bs; /* block col. index */ 2066 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2067 } 2068 rowcount++; 2069 } 2070 s_browlengths[i] += nmask; 2071 2072 /* zero out the mask elements we set */ 2073 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2074 } 2075 2076 /* create our matrix */ 2077 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2078 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2079 ierr = MatSetType(B,type);CHKERRQ(ierr); 2080 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2081 a = (Mat_SeqSBAIJ*)B->data; 2082 2083 /* set matrix "i" values */ 2084 a->i[0] = 0; 2085 for (i=1; i<= mbs; i++) { 2086 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2087 a->ilen[i-1] = s_browlengths[i-1]; 2088 } 2089 a->nz = a->i[mbs]; 2090 2091 /* read in nonzero values */ 2092 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2093 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2094 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2095 2096 /* set "a" and "j" values into matrix */ 2097 nzcount = 0; jcount = 0; 2098 for (i=0; i<mbs; i++) { 2099 nzcountb = nzcount; 2100 nmask = 0; 2101 for (j=0; j<bs; j++) { 2102 kmax = rowlengths[i*bs+j]; 2103 for (k=0; k<kmax; k++) { 2104 tmp = jj[nzcount++]/bs; /* block col. index */ 2105 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2106 } 2107 } 2108 /* sort the masked values */ 2109 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2110 2111 /* set "j" values into matrix */ 2112 maskcount = 1; 2113 for (j=0; j<nmask; j++) { 2114 a->j[jcount++] = masked[j]; 2115 mask[masked[j]] = maskcount++; 2116 } 2117 2118 /* set "a" values into matrix */ 2119 ishift = bs2*a->i[i]; 2120 for (j=0; j<bs; j++) { 2121 kmax = rowlengths[i*bs+j]; 2122 for (k=0; k<kmax; k++) { 2123 tmp = jj[nzcountb]/bs ; /* block col. index */ 2124 if (tmp >= i){ 2125 block = mask[tmp] - 1; 2126 point = jj[nzcountb] - bs*tmp; 2127 idx = ishift + bs2*block + j + bs*point; 2128 a->a[idx] = aa[nzcountb]; 2129 } 2130 nzcountb++; 2131 } 2132 } 2133 /* zero out the mask elements we set */ 2134 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2135 } 2136 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2137 2138 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2139 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2140 ierr = PetscFree(aa);CHKERRQ(ierr); 2141 ierr = PetscFree(jj);CHKERRQ(ierr); 2142 ierr = PetscFree(mask);CHKERRQ(ierr); 2143 2144 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2145 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2146 ierr = MatView_Private(B);CHKERRQ(ierr); 2147 *A = B; 2148 PetscFunctionReturn(0); 2149 } 2150 2151 2152 2153 #undef __FUNCT__ 2154 #define __FUNCT__ "MatRelax_SeqSBAIJ" 2155 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 2156 { 2157 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 2158 const MatScalar *aa=a->a,*v,*v1; 2159 PetscScalar *x,*b,*t,sum,d; 2160 MatScalar tmp; 2161 PetscErrorCode ierr; 2162 PetscInt m=a->mbs,bs=A->rmap->bs,j; 2163 const PetscInt *ai=a->i,*aj=a->j,*vj,*vj1; 2164 PetscInt nz,nz1,i; 2165 2166 PetscFunctionBegin; 2167 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 2168 2169 its = its*lits; 2170 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 2171 2172 if (bs > 1) SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 2173 2174 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 2175 if (xx != bb) { 2176 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2177 } else { 2178 b = x; 2179 } 2180 2181 if (!a->idiagvalid) { 2182 if (!a->idiag) { 2183 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->idiag);CHKERRQ(ierr); 2184 } 2185 for (i=0; i<a->mbs; i++) a->idiag[i] = 1.0/a->a[a->i[i]]; 2186 a->idiagvalid = PETSC_TRUE; 2187 } 2188 2189 if (!a->relax_work) { 2190 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2191 } 2192 t = a->relax_work; 2193 2194 2195 if (flag & SOR_ZERO_INITIAL_GUESS) { 2196 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2197 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2198 2199 for (i=0; i<m; i++){ 2200 v = aa + ai[i] + 1; 2201 vj = aj + ai[i] + 1; 2202 nz = ai[i+1] - ai[i] - 1; 2203 x[i] = tmp = omega*t[i]*a->idiag[i]; 2204 for (j=0; j<nz; j++) { 2205 t[vj[j]] -= tmp*v[j]; 2206 } 2207 } 2208 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2209 } 2210 2211 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2212 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2213 t = b; 2214 } 2215 2216 for (i=m-1; i>=0; i--){ 2217 v = aa + ai[i] + 1; 2218 vj = aj + ai[i] + 1; 2219 nz = ai[i+1] - ai[i] - 1; 2220 sum = t[i]; 2221 PetscSparseDenseMinusDot(sum,x,v,vj,nz); 2222 x[i] = (1-omega)*x[i] + omega*sum*a->idiag[i]; 2223 } 2224 t = a->relax_work; 2225 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 2226 } 2227 its--; 2228 } 2229 2230 while (its--) { 2231 /* 2232 forward sweep: 2233 for i=0,...,m-1: 2234 sum[i] = (b[i] - U(i,:)x )/d[i]; 2235 x[i] = (1-omega)x[i] + omega*sum[i]; 2236 b = b - x[i]*U^T(i,:); 2237 2238 */ 2239 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2240 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2241 2242 for (i=0; i<m; i++){ 2243 d = *(aa + ai[i]); /* diag[i] */ 2244 v = aa + ai[i] + 1; v1=v; 2245 vj = aj + ai[i] + 1; vj1=vj; 2246 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2247 sum = t[i]; 2248 ierr = PetscLogFlops(4.0*nz-2);CHKERRQ(ierr); 2249 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2250 x[i] = (1-omega)*x[i] + omega*sum/d; 2251 while (nz--) t[*vj++] -= x[i]*(*v++); 2252 } 2253 } 2254 2255 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2256 /* 2257 backward sweep: 2258 b = b - x[i]*U^T(i,:), i=0,...,n-2 2259 for i=m-1,...,0: 2260 sum[i] = (b[i] - U(i,:)x )/d[i]; 2261 x[i] = (1-omega)x[i] + omega*sum[i]; 2262 */ 2263 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2264 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2265 2266 for (i=0; i<m-1; i++){ /* update rhs */ 2267 v = aa + ai[i] + 1; 2268 vj = aj + ai[i] + 1; 2269 nz = ai[i+1] - ai[i] - 1; 2270 ierr = PetscLogFlops(2.0*nz-1);CHKERRQ(ierr); 2271 while (nz--) t[*vj++] -= x[i]*(*v++); 2272 } 2273 for (i=m-1; i>=0; i--){ 2274 d = *(aa + ai[i]); 2275 v = aa + ai[i] + 1; 2276 vj = aj + ai[i] + 1; 2277 nz = ai[i+1] - ai[i] - 1; 2278 ierr = PetscLogFlops(2.0*nz-1);CHKERRQ(ierr); 2279 sum = t[i]; 2280 while (nz--) sum -= x[*vj++]*(*v++); 2281 x[i] = (1-omega)*x[i] + omega*sum/d; 2282 } 2283 } 2284 } 2285 2286 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2287 if (bb != xx) { 2288 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2289 } 2290 PetscFunctionReturn(0); 2291 } 2292 2293 #undef __FUNCT__ 2294 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2295 /*@ 2296 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2297 (upper triangular entries in CSR format) provided by the user. 2298 2299 Collective on MPI_Comm 2300 2301 Input Parameters: 2302 + comm - must be an MPI communicator of size 1 2303 . bs - size of block 2304 . m - number of rows 2305 . n - number of columns 2306 . i - row indices 2307 . j - column indices 2308 - a - matrix values 2309 2310 Output Parameter: 2311 . mat - the matrix 2312 2313 Level: intermediate 2314 2315 Notes: 2316 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2317 once the matrix is destroyed 2318 2319 You cannot set new nonzero locations into this matrix, that will generate an error. 2320 2321 The i and j indices are 0 based 2322 2323 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2324 2325 @*/ 2326 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2327 { 2328 PetscErrorCode ierr; 2329 PetscInt ii; 2330 Mat_SeqSBAIJ *sbaij; 2331 2332 PetscFunctionBegin; 2333 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2334 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2335 2336 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2337 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2338 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2339 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2340 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2341 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2342 2343 sbaij->i = i; 2344 sbaij->j = j; 2345 sbaij->a = a; 2346 sbaij->singlemalloc = PETSC_FALSE; 2347 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2348 sbaij->free_a = PETSC_FALSE; 2349 sbaij->free_ij = PETSC_FALSE; 2350 2351 for (ii=0; ii<m; ii++) { 2352 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2353 #if defined(PETSC_USE_DEBUG) 2354 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]); 2355 #endif 2356 } 2357 #if defined(PETSC_USE_DEBUG) 2358 for (ii=0; ii<sbaij->i[m]; ii++) { 2359 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2360 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2361 } 2362 #endif 2363 2364 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2365 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2366 PetscFunctionReturn(0); 2367 } 2368 2369 2370 2371 2372 2373