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 ierr = MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew);CHKERRQ(ierr); 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 ierr = MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew);CHKERRQ(ierr); 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 #if defined(PETSC_HAVE_CHOLMOD) 1734 extern PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1735 #endif 1736 EXTERN_C_END 1737 1738 /*MC 1739 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1740 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1741 1742 Options Database Keys: 1743 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1744 1745 Level: beginner 1746 1747 .seealso: MatCreateSeqSBAIJ 1748 M*/ 1749 1750 EXTERN_C_BEGIN 1751 #undef __FUNCT__ 1752 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1753 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1754 { 1755 Mat_SeqSBAIJ *b; 1756 PetscErrorCode ierr; 1757 PetscMPIInt size; 1758 PetscTruth no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1759 1760 PetscFunctionBegin; 1761 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1762 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1763 1764 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1765 B->data = (void*)b; 1766 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1767 B->ops->destroy = MatDestroy_SeqSBAIJ; 1768 B->ops->view = MatView_SeqSBAIJ; 1769 B->mapping = 0; 1770 b->row = 0; 1771 b->icol = 0; 1772 b->reallocs = 0; 1773 b->saved_values = 0; 1774 b->inode.limit = 5; 1775 b->inode.max_limit = 5; 1776 1777 b->roworiented = PETSC_TRUE; 1778 b->nonew = 0; 1779 b->diag = 0; 1780 b->solve_work = 0; 1781 b->mult_work = 0; 1782 B->spptr = 0; 1783 B->info.nz_unneeded = (PetscReal)b->maxnz; 1784 b->keepnonzeropattern = PETSC_FALSE; 1785 b->xtoy = 0; 1786 b->XtoY = 0; 1787 1788 b->inew = 0; 1789 b->jnew = 0; 1790 b->anew = 0; 1791 b->a2anew = 0; 1792 b->permute = PETSC_FALSE; 1793 1794 b->ignore_ltriangular = PETSC_FALSE; 1795 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,PETSC_NULL);CHKERRQ(ierr); 1796 1797 b->getrow_utriangular = PETSC_FALSE; 1798 ierr = PetscOptionsGetTruth(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,PETSC_NULL);CHKERRQ(ierr); 1799 1800 #if defined(PETSC_HAVE_PASTIX) 1801 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C", 1802 "MatGetFactor_seqsbaij_pastix", 1803 MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1804 #endif 1805 #if defined(PETSC_HAVE_SPOOLES) 1806 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_spooles_C", 1807 "MatGetFactor_seqsbaij_spooles", 1808 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1809 #endif 1810 #if defined(PETSC_HAVE_MUMPS) 1811 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C", 1812 "MatGetFactor_sbaij_mumps", 1813 MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1814 #endif 1815 #if defined(PETSC_HAVE_CHOLMOD) 1816 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C", 1817 "MatGetFactor_seqsbaij_cholmod", 1818 MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1819 #endif 1820 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C", 1821 "MatGetFactorAvailable_seqsbaij_petsc", 1822 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1823 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C", 1824 "MatGetFactor_seqsbaij_petsc", 1825 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1826 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1827 "MatStoreValues_SeqSBAIJ", 1828 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1829 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1830 "MatRetrieveValues_SeqSBAIJ", 1831 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1832 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1833 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1834 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1835 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1836 "MatConvert_SeqSBAIJ_SeqAIJ", 1837 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1838 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1839 "MatConvert_SeqSBAIJ_SeqBAIJ", 1840 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1841 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1842 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1843 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1844 1845 B->symmetric = PETSC_TRUE; 1846 B->structurally_symmetric = PETSC_TRUE; 1847 B->symmetric_set = PETSC_TRUE; 1848 B->structurally_symmetric_set = PETSC_TRUE; 1849 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1850 1851 ierr = PetscOptionsBegin(((PetscObject)B)->comm,((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1852 ierr = PetscOptionsTruth("-mat_no_unroll","Do not optimize for inodes (slower)",PETSC_NULL,no_unroll,&no_unroll,PETSC_NULL);CHKERRQ(ierr); 1853 if (no_unroll) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr);} 1854 ierr = PetscOptionsTruth("-mat_no_inode","Do not optimize for inodes (slower)",PETSC_NULL,no_inode,&no_inode,PETSC_NULL);CHKERRQ(ierr); 1855 if (no_inode) {ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr);} 1856 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); 1857 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1858 b->inode.use = (PetscTruth)(!(no_unroll || no_inode)); 1859 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1860 1861 PetscFunctionReturn(0); 1862 } 1863 EXTERN_C_END 1864 1865 #undef __FUNCT__ 1866 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1867 /*@C 1868 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1869 compressed row) format. For good matrix assembly performance the 1870 user should preallocate the matrix storage by setting the parameter nz 1871 (or the array nnz). By setting these parameters accurately, performance 1872 during matrix assembly can be increased by more than a factor of 50. 1873 1874 Collective on Mat 1875 1876 Input Parameters: 1877 + A - the symmetric matrix 1878 . bs - size of block 1879 . nz - number of block nonzeros per block row (same for all rows) 1880 - nnz - array containing the number of block nonzeros in the upper triangular plus 1881 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1882 1883 Options Database Keys: 1884 . -mat_no_unroll - uses code that does not unroll the loops in the 1885 block calculations (much slower) 1886 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1887 1888 Level: intermediate 1889 1890 Notes: 1891 Specify the preallocated storage with either nz or nnz (not both). 1892 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1893 allocation. For additional details, see the users manual chapter on 1894 matrices. 1895 1896 You can call MatGetInfo() to get information on how effective the preallocation was; 1897 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1898 You can also run with the option -info and look for messages with the string 1899 malloc in them to see if additional memory allocation was needed. 1900 1901 If the nnz parameter is given then the nz parameter is ignored 1902 1903 1904 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1905 @*/ 1906 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1907 { 1908 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1909 1910 PetscFunctionBegin; 1911 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1912 if (f) { 1913 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1914 } 1915 PetscFunctionReturn(0); 1916 } 1917 1918 #undef __FUNCT__ 1919 #define __FUNCT__ "MatCreateSeqSBAIJ" 1920 /*@C 1921 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1922 compressed row) format. For good matrix assembly performance the 1923 user should preallocate the matrix storage by setting the parameter nz 1924 (or the array nnz). By setting these parameters accurately, performance 1925 during matrix assembly can be increased by more than a factor of 50. 1926 1927 Collective on MPI_Comm 1928 1929 Input Parameters: 1930 + comm - MPI communicator, set to PETSC_COMM_SELF 1931 . bs - size of block 1932 . m - number of rows, or number of columns 1933 . nz - number of block nonzeros per block row (same for all rows) 1934 - nnz - array containing the number of block nonzeros in the upper triangular plus 1935 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1936 1937 Output Parameter: 1938 . A - the symmetric matrix 1939 1940 Options Database Keys: 1941 . -mat_no_unroll - uses code that does not unroll the loops in the 1942 block calculations (much slower) 1943 . -mat_block_size - size of the blocks to use 1944 1945 Level: intermediate 1946 1947 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1948 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1949 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1950 1951 Notes: 1952 The number of rows and columns must be divisible by blocksize. 1953 This matrix type does not support complex Hermitian operation. 1954 1955 Specify the preallocated storage with either nz or nnz (not both). 1956 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1957 allocation. For additional details, see the users manual chapter on 1958 matrices. 1959 1960 If the nnz parameter is given then the nz parameter is ignored 1961 1962 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1963 @*/ 1964 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1965 { 1966 PetscErrorCode ierr; 1967 1968 PetscFunctionBegin; 1969 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1970 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1971 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1972 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1973 PetscFunctionReturn(0); 1974 } 1975 1976 #undef __FUNCT__ 1977 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1978 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1979 { 1980 Mat C; 1981 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1982 PetscErrorCode ierr; 1983 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1984 1985 PetscFunctionBegin; 1986 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 1987 1988 *B = 0; 1989 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1990 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1991 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 1992 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1993 c = (Mat_SeqSBAIJ*)C->data; 1994 1995 C->preallocated = PETSC_TRUE; 1996 C->factortype = A->factortype; 1997 c->row = 0; 1998 c->icol = 0; 1999 c->saved_values = 0; 2000 c->keepnonzeropattern = a->keepnonzeropattern; 2001 C->assembled = PETSC_TRUE; 2002 2003 ierr = PetscLayoutCopy(A->rmap,&C->rmap);CHKERRQ(ierr); 2004 ierr = PetscLayoutCopy(A->cmap,&C->cmap);CHKERRQ(ierr); 2005 c->bs2 = a->bs2; 2006 c->mbs = a->mbs; 2007 c->nbs = a->nbs; 2008 2009 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2010 c->imax = a->imax; 2011 c->ilen = a->ilen; 2012 c->free_imax_ilen = PETSC_FALSE; 2013 } else { 2014 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2015 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2016 for (i=0; i<mbs; i++) { 2017 c->imax[i] = a->imax[i]; 2018 c->ilen[i] = a->ilen[i]; 2019 } 2020 c->free_imax_ilen = PETSC_TRUE; 2021 } 2022 2023 /* allocate the matrix space */ 2024 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2025 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2026 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2027 c->singlemalloc = PETSC_FALSE; 2028 c->free_ij = PETSC_FALSE; 2029 c->parent = A; 2030 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2031 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2032 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2033 } else { 2034 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2035 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2036 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2037 c->singlemalloc = PETSC_TRUE; 2038 c->free_ij = PETSC_TRUE; 2039 } 2040 if (mbs > 0) { 2041 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2042 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2043 } 2044 if (cpvalues == MAT_COPY_VALUES) { 2045 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2046 } else { 2047 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2048 } 2049 if (a->jshort) { 2050 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2051 c->jshort = a->jshort; 2052 c->free_jshort = PETSC_FALSE; 2053 } else { 2054 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2055 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2056 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2057 c->free_jshort = PETSC_TRUE; 2058 } 2059 } 2060 } 2061 2062 c->roworiented = a->roworiented; 2063 c->nonew = a->nonew; 2064 2065 if (a->diag) { 2066 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2067 c->diag = a->diag; 2068 c->free_diag = PETSC_FALSE; 2069 } else { 2070 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2071 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2072 for (i=0; i<mbs; i++) { 2073 c->diag[i] = a->diag[i]; 2074 } 2075 c->free_diag = PETSC_TRUE; 2076 } 2077 } else c->diag = 0; 2078 c->nz = a->nz; 2079 c->maxnz = bs2*a->nz; /* Since we allocate exactly the right amount */ 2080 c->solve_work = 0; 2081 c->mult_work = 0; 2082 c->free_a = PETSC_TRUE; 2083 *B = C; 2084 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2085 PetscFunctionReturn(0); 2086 } 2087 2088 #undef __FUNCT__ 2089 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2090 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 2091 { 2092 Mat_SeqSBAIJ *a; 2093 Mat B; 2094 PetscErrorCode ierr; 2095 int fd; 2096 PetscMPIInt size; 2097 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2098 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2099 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 2100 PetscInt *masked,nmask,tmp,bs2,ishift; 2101 PetscScalar *aa; 2102 MPI_Comm comm = ((PetscObject)viewer)->comm; 2103 2104 PetscFunctionBegin; 2105 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2106 bs2 = bs*bs; 2107 2108 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2109 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2110 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2111 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2112 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2113 M = header[1]; N = header[2]; nz = header[3]; 2114 2115 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2116 2117 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2118 2119 /* 2120 This code adds extra rows to make sure the number of rows is 2121 divisible by the blocksize 2122 */ 2123 mbs = M/bs; 2124 extra_rows = bs - M + bs*(mbs); 2125 if (extra_rows == bs) extra_rows = 0; 2126 else mbs++; 2127 if (extra_rows) { 2128 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2129 } 2130 2131 /* read in row lengths */ 2132 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2133 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2134 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2135 2136 /* read in column indices */ 2137 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2138 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2139 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2140 2141 /* loop over row lengths determining block row lengths */ 2142 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2143 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2144 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2145 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2146 rowcount = 0; 2147 nzcount = 0; 2148 for (i=0; i<mbs; i++) { 2149 nmask = 0; 2150 for (j=0; j<bs; j++) { 2151 kmax = rowlengths[rowcount]; 2152 for (k=0; k<kmax; k++) { 2153 tmp = jj[nzcount++]/bs; /* block col. index */ 2154 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2155 } 2156 rowcount++; 2157 } 2158 s_browlengths[i] += nmask; 2159 2160 /* zero out the mask elements we set */ 2161 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2162 } 2163 2164 /* create our matrix */ 2165 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 2166 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2167 ierr = MatSetType(B,type);CHKERRQ(ierr); 2168 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 2169 a = (Mat_SeqSBAIJ*)B->data; 2170 2171 /* set matrix "i" values */ 2172 a->i[0] = 0; 2173 for (i=1; i<= mbs; i++) { 2174 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2175 a->ilen[i-1] = s_browlengths[i-1]; 2176 } 2177 a->nz = a->i[mbs]; 2178 2179 /* read in nonzero values */ 2180 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2181 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2182 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2183 2184 /* set "a" and "j" values into matrix */ 2185 nzcount = 0; jcount = 0; 2186 for (i=0; i<mbs; i++) { 2187 nzcountb = nzcount; 2188 nmask = 0; 2189 for (j=0; j<bs; j++) { 2190 kmax = rowlengths[i*bs+j]; 2191 for (k=0; k<kmax; k++) { 2192 tmp = jj[nzcount++]/bs; /* block col. index */ 2193 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2194 } 2195 } 2196 /* sort the masked values */ 2197 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2198 2199 /* set "j" values into matrix */ 2200 maskcount = 1; 2201 for (j=0; j<nmask; j++) { 2202 a->j[jcount++] = masked[j]; 2203 mask[masked[j]] = maskcount++; 2204 } 2205 2206 /* set "a" values into matrix */ 2207 ishift = bs2*a->i[i]; 2208 for (j=0; j<bs; j++) { 2209 kmax = rowlengths[i*bs+j]; 2210 for (k=0; k<kmax; k++) { 2211 tmp = jj[nzcountb]/bs ; /* block col. index */ 2212 if (tmp >= i){ 2213 block = mask[tmp] - 1; 2214 point = jj[nzcountb] - bs*tmp; 2215 idx = ishift + bs2*block + j + bs*point; 2216 a->a[idx] = aa[nzcountb]; 2217 } 2218 nzcountb++; 2219 } 2220 } 2221 /* zero out the mask elements we set */ 2222 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2223 } 2224 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2225 2226 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2227 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2228 ierr = PetscFree(aa);CHKERRQ(ierr); 2229 ierr = PetscFree(jj);CHKERRQ(ierr); 2230 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2231 2232 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2233 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2234 ierr = MatView_Private(B);CHKERRQ(ierr); 2235 *A = B; 2236 PetscFunctionReturn(0); 2237 } 2238 2239 #undef __FUNCT__ 2240 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2241 /*@ 2242 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2243 (upper triangular entries in CSR format) provided by the user. 2244 2245 Collective on MPI_Comm 2246 2247 Input Parameters: 2248 + comm - must be an MPI communicator of size 1 2249 . bs - size of block 2250 . m - number of rows 2251 . n - number of columns 2252 . i - row indices 2253 . j - column indices 2254 - a - matrix values 2255 2256 Output Parameter: 2257 . mat - the matrix 2258 2259 Level: intermediate 2260 2261 Notes: 2262 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2263 once the matrix is destroyed 2264 2265 You cannot set new nonzero locations into this matrix, that will generate an error. 2266 2267 The i and j indices are 0 based 2268 2269 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2270 2271 @*/ 2272 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2273 { 2274 PetscErrorCode ierr; 2275 PetscInt ii; 2276 Mat_SeqSBAIJ *sbaij; 2277 2278 PetscFunctionBegin; 2279 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2280 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2281 2282 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2283 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2284 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2285 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2286 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2287 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2288 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2289 2290 sbaij->i = i; 2291 sbaij->j = j; 2292 sbaij->a = a; 2293 sbaij->singlemalloc = PETSC_FALSE; 2294 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2295 sbaij->free_a = PETSC_FALSE; 2296 sbaij->free_ij = PETSC_FALSE; 2297 2298 for (ii=0; ii<m; ii++) { 2299 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2300 #if defined(PETSC_USE_DEBUG) 2301 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]); 2302 #endif 2303 } 2304 #if defined(PETSC_USE_DEBUG) 2305 for (ii=0; ii<sbaij->i[m]; ii++) { 2306 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2307 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]); 2308 } 2309 #endif 2310 2311 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2312 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2313 PetscFunctionReturn(0); 2314 } 2315 2316 2317 2318 2319 2320