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