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