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