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