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