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