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