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