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