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*bs2; 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; 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. For additional details, see the users manual chapter on 1903 matrices. 1904 1905 You can call MatGetInfo() to get information on how effective the preallocation was; 1906 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1907 You can also run with the option -info and look for messages with the string 1908 malloc in them to see if additional memory allocation was needed. 1909 1910 If the nnz parameter is given then the nz parameter is ignored 1911 1912 1913 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1914 @*/ 1915 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1916 { 1917 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1918 1919 PetscFunctionBegin; 1920 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1921 if (f) { 1922 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1923 } 1924 PetscFunctionReturn(0); 1925 } 1926 1927 #undef __FUNCT__ 1928 #define __FUNCT__ "MatCreateSeqSBAIJ" 1929 /*@C 1930 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1931 compressed row) format. For good matrix assembly performance the 1932 user should preallocate the matrix storage by setting the parameter nz 1933 (or the array nnz). By setting these parameters accurately, performance 1934 during matrix assembly can be increased by more than a factor of 50. 1935 1936 Collective on MPI_Comm 1937 1938 Input Parameters: 1939 + comm - MPI communicator, set to PETSC_COMM_SELF 1940 . bs - size of block 1941 . m - number of rows, or number of columns 1942 . nz - number of block nonzeros per block row (same for all rows) 1943 - nnz - array containing the number of block nonzeros in the upper triangular plus 1944 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1945 1946 Output Parameter: 1947 . A - the symmetric matrix 1948 1949 Options Database Keys: 1950 . -mat_no_unroll - uses code that does not unroll the loops in the 1951 block calculations (much slower) 1952 . -mat_block_size - size of the blocks to use 1953 1954 Level: intermediate 1955 1956 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1957 MatXXXXSetPreallocation() paradgm instead of this routine directly. 1958 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1959 1960 Notes: 1961 The number of rows and columns must be divisible by blocksize. 1962 This matrix type does not support complex Hermitian operation. 1963 1964 Specify the preallocated storage with either nz or nnz (not both). 1965 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1966 allocation. For additional details, see the users manual chapter on 1967 matrices. 1968 1969 If the nnz parameter is given then the nz parameter is ignored 1970 1971 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1972 @*/ 1973 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1974 { 1975 PetscErrorCode ierr; 1976 1977 PetscFunctionBegin; 1978 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1979 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1980 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1981 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1982 PetscFunctionReturn(0); 1983 } 1984 1985 #undef __FUNCT__ 1986 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1987 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1988 { 1989 Mat C; 1990 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1991 PetscErrorCode ierr; 1992 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1993 1994 PetscFunctionBegin; 1995 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 1996 1997 *B = 0; 1998 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1999 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2000 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2001 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2002 c = (Mat_SeqSBAIJ*)C->data; 2003 2004 C->preallocated = PETSC_TRUE; 2005 C->factortype = A->factortype; 2006 c->row = 0; 2007 c->icol = 0; 2008 c->saved_values = 0; 2009 c->keepnonzeropattern = a->keepnonzeropattern; 2010 C->assembled = PETSC_TRUE; 2011 2012 ierr = PetscLayoutCopy(A->rmap,&C->rmap);CHKERRQ(ierr); 2013 ierr = PetscLayoutCopy(A->cmap,&C->cmap);CHKERRQ(ierr); 2014 c->bs2 = a->bs2; 2015 c->mbs = a->mbs; 2016 c->nbs = a->nbs; 2017 2018 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2019 c->imax = a->imax; 2020 c->ilen = a->ilen; 2021 c->free_imax_ilen = PETSC_FALSE; 2022 } else { 2023 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2024 ierr = PetscLogObjectMemory(C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2025 for (i=0; i<mbs; i++) { 2026 c->imax[i] = a->imax[i]; 2027 c->ilen[i] = a->ilen[i]; 2028 } 2029 c->free_imax_ilen = PETSC_TRUE; 2030 } 2031 2032 /* allocate the matrix space */ 2033 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2034 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2035 ierr = PetscLogObjectMemory(C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2036 c->singlemalloc = PETSC_FALSE; 2037 c->free_ij = PETSC_FALSE; 2038 c->parent = A; 2039 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2040 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2041 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2042 } else { 2043 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2044 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2045 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2046 c->singlemalloc = PETSC_TRUE; 2047 c->free_ij = PETSC_TRUE; 2048 } 2049 if (mbs > 0) { 2050 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2051 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2052 } 2053 if (cpvalues == MAT_COPY_VALUES) { 2054 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2055 } else { 2056 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2057 } 2058 if (a->jshort) { 2059 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2060 c->jshort = a->jshort; 2061 c->free_jshort = PETSC_FALSE; 2062 } else { 2063 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2064 ierr = PetscLogObjectMemory(C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2065 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2066 c->free_jshort = PETSC_TRUE; 2067 } 2068 } 2069 } 2070 2071 c->roworiented = a->roworiented; 2072 c->nonew = a->nonew; 2073 2074 if (a->diag) { 2075 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2076 c->diag = a->diag; 2077 c->free_diag = PETSC_FALSE; 2078 } else { 2079 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2080 ierr = PetscLogObjectMemory(C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2081 for (i=0; i<mbs; i++) { 2082 c->diag[i] = a->diag[i]; 2083 } 2084 c->free_diag = PETSC_TRUE; 2085 } 2086 } else c->diag = 0; 2087 c->nz = a->nz; 2088 c->maxnz = bs2*a->nz; /* Since we allocate exactly the right amount */ 2089 c->solve_work = 0; 2090 c->mult_work = 0; 2091 c->free_a = PETSC_TRUE; 2092 *B = C; 2093 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2094 PetscFunctionReturn(0); 2095 } 2096 2097 #undef __FUNCT__ 2098 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2099 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, Mat newmat) 2100 { 2101 Mat_SeqSBAIJ *a; 2102 PetscErrorCode ierr; 2103 int fd; 2104 PetscMPIInt size; 2105 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2106 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2107 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2108 PetscInt *masked,nmask,tmp,bs2,ishift; 2109 PetscScalar *aa; 2110 MPI_Comm comm = ((PetscObject)viewer)->comm; 2111 2112 PetscFunctionBegin; 2113 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 2114 bs2 = bs*bs; 2115 2116 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2117 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2118 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2119 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2120 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2121 M = header[1]; N = header[2]; nz = header[3]; 2122 2123 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2124 2125 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2126 2127 /* 2128 This code adds extra rows to make sure the number of rows is 2129 divisible by the blocksize 2130 */ 2131 mbs = M/bs; 2132 extra_rows = bs - M + bs*(mbs); 2133 if (extra_rows == bs) extra_rows = 0; 2134 else mbs++; 2135 if (extra_rows) { 2136 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2137 } 2138 2139 /* Set global sizes if not already set */ 2140 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2141 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2142 } else { /* Check if the matrix global sizes are correct */ 2143 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2144 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); 2145 } 2146 2147 /* read in row lengths */ 2148 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2149 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2150 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2151 2152 /* read in column indices */ 2153 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2154 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2155 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2156 2157 /* loop over row lengths determining block row lengths */ 2158 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2159 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2160 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2161 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2162 rowcount = 0; 2163 nzcount = 0; 2164 for (i=0; i<mbs; i++) { 2165 nmask = 0; 2166 for (j=0; j<bs; j++) { 2167 kmax = rowlengths[rowcount]; 2168 for (k=0; k<kmax; k++) { 2169 tmp = jj[nzcount++]/bs; /* block col. index */ 2170 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2171 } 2172 rowcount++; 2173 } 2174 s_browlengths[i] += nmask; 2175 2176 /* zero out the mask elements we set */ 2177 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2178 } 2179 2180 /* Do preallocation */ 2181 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2182 a = (Mat_SeqSBAIJ*)newmat->data; 2183 2184 /* set matrix "i" values */ 2185 a->i[0] = 0; 2186 for (i=1; i<= mbs; i++) { 2187 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2188 a->ilen[i-1] = s_browlengths[i-1]; 2189 } 2190 a->nz = a->i[mbs]; 2191 2192 /* read in nonzero values */ 2193 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2194 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2195 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2196 2197 /* set "a" and "j" values into matrix */ 2198 nzcount = 0; jcount = 0; 2199 for (i=0; i<mbs; i++) { 2200 nzcountb = nzcount; 2201 nmask = 0; 2202 for (j=0; j<bs; j++) { 2203 kmax = rowlengths[i*bs+j]; 2204 for (k=0; k<kmax; k++) { 2205 tmp = jj[nzcount++]/bs; /* block col. index */ 2206 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2207 } 2208 } 2209 /* sort the masked values */ 2210 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2211 2212 /* set "j" values into matrix */ 2213 maskcount = 1; 2214 for (j=0; j<nmask; j++) { 2215 a->j[jcount++] = masked[j]; 2216 mask[masked[j]] = maskcount++; 2217 } 2218 2219 /* set "a" values into matrix */ 2220 ishift = bs2*a->i[i]; 2221 for (j=0; j<bs; j++) { 2222 kmax = rowlengths[i*bs+j]; 2223 for (k=0; k<kmax; k++) { 2224 tmp = jj[nzcountb]/bs ; /* block col. index */ 2225 if (tmp >= i){ 2226 block = mask[tmp] - 1; 2227 point = jj[nzcountb] - bs*tmp; 2228 idx = ishift + bs2*block + j + bs*point; 2229 a->a[idx] = aa[nzcountb]; 2230 } 2231 nzcountb++; 2232 } 2233 } 2234 /* zero out the mask elements we set */ 2235 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2236 } 2237 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2238 2239 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2240 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2241 ierr = PetscFree(aa);CHKERRQ(ierr); 2242 ierr = PetscFree(jj);CHKERRQ(ierr); 2243 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2244 2245 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2246 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2247 ierr = MatView_Private(newmat);CHKERRQ(ierr); 2248 PetscFunctionReturn(0); 2249 } 2250 2251 #undef __FUNCT__ 2252 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2253 /*@ 2254 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2255 (upper triangular entries in CSR format) provided by the user. 2256 2257 Collective on MPI_Comm 2258 2259 Input Parameters: 2260 + comm - must be an MPI communicator of size 1 2261 . bs - size of block 2262 . m - number of rows 2263 . n - number of columns 2264 . i - row indices 2265 . j - column indices 2266 - a - matrix values 2267 2268 Output Parameter: 2269 . mat - the matrix 2270 2271 Level: intermediate 2272 2273 Notes: 2274 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2275 once the matrix is destroyed 2276 2277 You cannot set new nonzero locations into this matrix, that will generate an error. 2278 2279 The i and j indices are 0 based 2280 2281 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2282 2283 @*/ 2284 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2285 { 2286 PetscErrorCode ierr; 2287 PetscInt ii; 2288 Mat_SeqSBAIJ *sbaij; 2289 2290 PetscFunctionBegin; 2291 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2292 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2293 2294 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2295 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2296 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2297 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2298 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2299 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2300 ierr = PetscLogObjectMemory(*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2301 2302 sbaij->i = i; 2303 sbaij->j = j; 2304 sbaij->a = a; 2305 sbaij->singlemalloc = PETSC_FALSE; 2306 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2307 sbaij->free_a = PETSC_FALSE; 2308 sbaij->free_ij = PETSC_FALSE; 2309 2310 for (ii=0; ii<m; ii++) { 2311 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2312 #if defined(PETSC_USE_DEBUG) 2313 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]); 2314 #endif 2315 } 2316 #if defined(PETSC_USE_DEBUG) 2317 for (ii=0; ii<sbaij->i[m]; ii++) { 2318 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2319 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]); 2320 } 2321 #endif 2322 2323 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2324 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2325 PetscFunctionReturn(0); 2326 } 2327 2328 2329 2330 2331 2332