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