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(inA,info,&outA);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 MatCholeskyFactor_SeqSBAIJ, 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 MatForwardSolve_SeqSBAIJ_N, 1182 MatBackwardSolve_SeqSBAIJ_N, 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_C_END 1547 1548 1549 EXTERN_C_BEGIN 1550 #undef __FUNCT__ 1551 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1552 PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1553 { 1554 PetscInt n = A->rmap->n; 1555 PetscErrorCode ierr; 1556 1557 PetscFunctionBegin; 1558 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 1559 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1560 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1561 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1562 ierr = MatSeqSBAIJSetPreallocation(*B,1,MAT_SKIP_ALLOCATION,PETSC_NULL);CHKERRQ(ierr); 1563 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1564 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1565 } else SETERRQ(PETSC_ERR_SUP,"Factor type not supported"); 1566 PetscFunctionReturn(0); 1567 } 1568 EXTERN_C_END 1569 1570 EXTERN_C_BEGIN 1571 #undef __FUNCT__ 1572 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1573 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscTruth *flg) 1574 { 1575 PetscFunctionBegin; 1576 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1577 *flg = PETSC_TRUE; 1578 } else { 1579 *flg = PETSC_FALSE; 1580 } 1581 PetscFunctionReturn(0); 1582 } 1583 EXTERN_C_END 1584 1585 EXTERN_C_BEGIN 1586 #if defined(PETSC_HAVE_MUMPS) 1587 extern PetscErrorCode MatGetFactor_seqsbaij_mumps(Mat,MatFactorType,Mat*); 1588 #endif 1589 #if defined(PETSC_HAVE_SPOOLES) 1590 extern PetscErrorCode MatGetFactor_seqsbaij_spooles(Mat,MatFactorType,Mat*); 1591 #endif 1592 EXTERN_C_END 1593 1594 /*MC 1595 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1596 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1597 1598 Options Database Keys: 1599 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1600 1601 Level: beginner 1602 1603 .seealso: MatCreateSeqSBAIJ 1604 M*/ 1605 1606 EXTERN_C_BEGIN 1607 #undef __FUNCT__ 1608 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1609 PetscErrorCode PETSCMAT_DLLEXPORT MatCreate_SeqSBAIJ(Mat B) 1610 { 1611 Mat_SeqSBAIJ *b; 1612 PetscErrorCode ierr; 1613 PetscMPIInt size; 1614 PetscTruth flg; 1615 1616 PetscFunctionBegin; 1617 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 1618 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1619 1620 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1621 B->data = (void*)b; 1622 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1623 B->ops->destroy = MatDestroy_SeqSBAIJ; 1624 B->ops->view = MatView_SeqSBAIJ; 1625 B->mapping = 0; 1626 b->row = 0; 1627 b->icol = 0; 1628 b->reallocs = 0; 1629 b->saved_values = 0; 1630 1631 1632 b->roworiented = PETSC_TRUE; 1633 b->nonew = 0; 1634 b->diag = 0; 1635 b->solve_work = 0; 1636 b->mult_work = 0; 1637 B->spptr = 0; 1638 b->keepzeroedrows = PETSC_FALSE; 1639 b->xtoy = 0; 1640 b->XtoY = 0; 1641 1642 b->inew = 0; 1643 b->jnew = 0; 1644 b->anew = 0; 1645 b->a2anew = 0; 1646 b->permute = PETSC_FALSE; 1647 1648 b->ignore_ltriangular = PETSC_FALSE; 1649 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_ignore_lower_triangular",&flg);CHKERRQ(ierr); 1650 if (flg) b->ignore_ltriangular = PETSC_TRUE; 1651 1652 b->getrow_utriangular = PETSC_FALSE; 1653 ierr = PetscOptionsHasName(PETSC_NULL,"-mat_getrow_uppertriangular",&flg);CHKERRQ(ierr); 1654 if (flg) b->getrow_utriangular = PETSC_TRUE; 1655 1656 #if defined(PETSC_HAVE_SPOOLES) 1657 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_spooles_C", 1658 "MatGetFactor_seqsbaij_spooles", 1659 MatGetFactor_seqsbaij_spooles);CHKERRQ(ierr); 1660 #endif 1661 #if defined(PETSC_HAVE_MUMPS) 1662 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_mumps_C", 1663 "MatGetFactor_seqsbaij_mumps", 1664 MatGetFactor_seqsbaij_mumps);CHKERRQ(ierr); 1665 #endif 1666 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_seqsbaij_petsc_C", 1667 "MatGetFactorAvailable_seqsbaij_petsc", 1668 MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1669 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_seqsbaij_petsc_C", 1670 "MatGetFactor_seqsbaij_petsc", 1671 MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1672 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C", 1673 "MatStoreValues_SeqSBAIJ", 1674 MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1675 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C", 1676 "MatRetrieveValues_SeqSBAIJ", 1677 (void*)MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1678 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C", 1679 "MatSeqSBAIJSetColumnIndices_SeqSBAIJ", 1680 MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1681 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqaij_C", 1682 "MatConvert_SeqSBAIJ_SeqAIJ", 1683 MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1684 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C", 1685 "MatConvert_SeqSBAIJ_SeqBAIJ", 1686 MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1687 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqSBAIJSetPreallocation_C", 1688 "MatSeqSBAIJSetPreallocation_SeqSBAIJ", 1689 MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1690 1691 B->symmetric = PETSC_TRUE; 1692 B->structurally_symmetric = PETSC_TRUE; 1693 B->symmetric_set = PETSC_TRUE; 1694 B->structurally_symmetric_set = PETSC_TRUE; 1695 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1696 PetscFunctionReturn(0); 1697 } 1698 EXTERN_C_END 1699 1700 #undef __FUNCT__ 1701 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 1702 /*@C 1703 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 1704 compressed row) format. For good matrix assembly performance the 1705 user should preallocate the matrix storage by setting the parameter nz 1706 (or the array nnz). By setting these parameters accurately, performance 1707 during matrix assembly can be increased by more than a factor of 50. 1708 1709 Collective on Mat 1710 1711 Input Parameters: 1712 + A - the symmetric matrix 1713 . bs - size of block 1714 . nz - number of block nonzeros per block row (same for all rows) 1715 - nnz - array containing the number of block nonzeros in the upper triangular plus 1716 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1717 1718 Options Database Keys: 1719 . -mat_no_unroll - uses code that does not unroll the loops in the 1720 block calculations (much slower) 1721 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 1722 1723 Level: intermediate 1724 1725 Notes: 1726 Specify the preallocated storage with either nz or nnz (not both). 1727 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1728 allocation. For additional details, see the users manual chapter on 1729 matrices. 1730 1731 You can call MatGetInfo() to get information on how effective the preallocation was; 1732 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 1733 You can also run with the option -info and look for messages with the string 1734 malloc in them to see if additional memory allocation was needed. 1735 1736 If the nnz parameter is given then the nz parameter is ignored 1737 1738 1739 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1740 @*/ 1741 PetscErrorCode PETSCMAT_DLLEXPORT MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 1742 { 1743 PetscErrorCode ierr,(*f)(Mat,PetscInt,PetscInt,const PetscInt[]); 1744 1745 PetscFunctionBegin; 1746 ierr = PetscObjectQueryFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",(void (**)(void))&f);CHKERRQ(ierr); 1747 if (f) { 1748 ierr = (*f)(B,bs,nz,nnz);CHKERRQ(ierr); 1749 } 1750 PetscFunctionReturn(0); 1751 } 1752 1753 #undef __FUNCT__ 1754 #define __FUNCT__ "MatCreateSeqSBAIJ" 1755 /*@C 1756 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 1757 compressed row) format. For good matrix assembly performance the 1758 user should preallocate the matrix storage by setting the parameter nz 1759 (or the array nnz). By setting these parameters accurately, performance 1760 during matrix assembly can be increased by more than a factor of 50. 1761 1762 Collective on MPI_Comm 1763 1764 Input Parameters: 1765 + comm - MPI communicator, set to PETSC_COMM_SELF 1766 . bs - size of block 1767 . m - number of rows, or number of columns 1768 . nz - number of block nonzeros per block row (same for all rows) 1769 - nnz - array containing the number of block nonzeros in the upper triangular plus 1770 diagonal portion of each block (possibly different for each block row) or PETSC_NULL 1771 1772 Output Parameter: 1773 . A - the symmetric matrix 1774 1775 Options Database Keys: 1776 . -mat_no_unroll - uses code that does not unroll the loops in the 1777 block calculations (much slower) 1778 . -mat_block_size - size of the blocks to use 1779 1780 Level: intermediate 1781 1782 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 1783 MatXXXXSetPreallocation() paradgm instead of this routine directly. This is definitely 1784 true if you plan to use the external direct solvers such as SuperLU, MUMPS or Spooles. 1785 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 1786 1787 Notes: 1788 The number of rows and columns must be divisible by blocksize. 1789 This matrix type does not support complex Hermitian operation. 1790 1791 Specify the preallocated storage with either nz or nnz (not both). 1792 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 1793 allocation. For additional details, see the users manual chapter on 1794 matrices. 1795 1796 If the nnz parameter is given then the nz parameter is ignored 1797 1798 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateMPISBAIJ() 1799 @*/ 1800 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 1801 { 1802 PetscErrorCode ierr; 1803 1804 PetscFunctionBegin; 1805 ierr = MatCreate(comm,A);CHKERRQ(ierr); 1806 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 1807 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 1808 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 1809 PetscFunctionReturn(0); 1810 } 1811 1812 #undef __FUNCT__ 1813 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 1814 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 1815 { 1816 Mat C; 1817 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 1818 PetscErrorCode ierr; 1819 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 1820 1821 PetscFunctionBegin; 1822 if (a->i[mbs] != nz) SETERRQ(PETSC_ERR_PLIB,"Corrupt matrix"); 1823 1824 *B = 0; 1825 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1826 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 1827 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1828 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 1829 c = (Mat_SeqSBAIJ*)C->data; 1830 1831 C->preallocated = PETSC_TRUE; 1832 C->factor = A->factor; 1833 c->row = 0; 1834 c->icol = 0; 1835 c->saved_values = 0; 1836 c->keepzeroedrows = a->keepzeroedrows; 1837 C->assembled = PETSC_TRUE; 1838 1839 ierr = PetscMapCopy(((PetscObject)A)->comm,A->rmap,C->rmap);CHKERRQ(ierr); 1840 ierr = PetscMapCopy(((PetscObject)A)->comm,A->cmap,C->cmap);CHKERRQ(ierr); 1841 c->bs2 = a->bs2; 1842 c->mbs = a->mbs; 1843 c->nbs = a->nbs; 1844 1845 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 1846 for (i=0; i<mbs; i++) { 1847 c->imax[i] = a->imax[i]; 1848 c->ilen[i] = a->ilen[i]; 1849 } 1850 1851 /* allocate the matrix space */ 1852 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 1853 c->singlemalloc = PETSC_TRUE; 1854 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1855 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 1856 if (mbs > 0) { 1857 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1858 if (cpvalues == MAT_COPY_VALUES) { 1859 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1860 } else { 1861 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 1862 } 1863 } 1864 1865 c->roworiented = a->roworiented; 1866 c->nonew = a->nonew; 1867 1868 if (a->diag) { 1869 ierr = PetscMalloc((mbs+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 1870 ierr = PetscLogObjectMemory(C,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 1871 for (i=0; i<mbs; i++) { 1872 c->diag[i] = a->diag[i]; 1873 } 1874 } else c->diag = 0; 1875 c->nz = a->nz; 1876 c->maxnz = a->maxnz; 1877 c->solve_work = 0; 1878 c->mult_work = 0; 1879 c->free_a = PETSC_TRUE; 1880 c->free_ij = PETSC_TRUE; 1881 *B = C; 1882 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 1883 PetscFunctionReturn(0); 1884 } 1885 1886 #undef __FUNCT__ 1887 #define __FUNCT__ "MatLoad_SeqSBAIJ" 1888 PetscErrorCode MatLoad_SeqSBAIJ(PetscViewer viewer, const MatType type,Mat *A) 1889 { 1890 Mat_SeqSBAIJ *a; 1891 Mat B; 1892 PetscErrorCode ierr; 1893 int fd; 1894 PetscMPIInt size; 1895 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 1896 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 1897 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows; 1898 PetscInt *masked,nmask,tmp,bs2,ishift; 1899 PetscScalar *aa; 1900 MPI_Comm comm = ((PetscObject)viewer)->comm; 1901 1902 PetscFunctionBegin; 1903 ierr = PetscOptionsGetInt(PETSC_NULL,"-matload_block_size",&bs,PETSC_NULL);CHKERRQ(ierr); 1904 bs2 = bs*bs; 1905 1906 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 1907 if (size > 1) SETERRQ(PETSC_ERR_ARG_WRONG,"view must have one processor"); 1908 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 1909 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 1910 if (header[0] != MAT_FILE_COOKIE) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 1911 M = header[1]; N = header[2]; nz = header[3]; 1912 1913 if (header[3] < 0) { 1914 SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 1915 } 1916 1917 if (M != N) SETERRQ(PETSC_ERR_SUP,"Can only do square matrices"); 1918 1919 /* 1920 This code adds extra rows to make sure the number of rows is 1921 divisible by the blocksize 1922 */ 1923 mbs = M/bs; 1924 extra_rows = bs - M + bs*(mbs); 1925 if (extra_rows == bs) extra_rows = 0; 1926 else mbs++; 1927 if (extra_rows) { 1928 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 1929 } 1930 1931 /* read in row lengths */ 1932 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 1933 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 1934 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 1935 1936 /* read in column indices */ 1937 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 1938 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 1939 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 1940 1941 /* loop over row lengths determining block row lengths */ 1942 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 1943 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1944 ierr = PetscMalloc(2*mbs*sizeof(PetscInt),&mask);CHKERRQ(ierr); 1945 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 1946 masked = mask + mbs; 1947 rowcount = 0; nzcount = 0; 1948 for (i=0; i<mbs; i++) { 1949 nmask = 0; 1950 for (j=0; j<bs; j++) { 1951 kmax = rowlengths[rowcount]; 1952 for (k=0; k<kmax; k++) { 1953 tmp = jj[nzcount++]/bs; /* block col. index */ 1954 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 1955 } 1956 rowcount++; 1957 } 1958 s_browlengths[i] += nmask; 1959 1960 /* zero out the mask elements we set */ 1961 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 1962 } 1963 1964 /* create our matrix */ 1965 ierr = MatCreate(comm,&B);CHKERRQ(ierr); 1966 ierr = MatSetSizes(B,M+extra_rows,N+extra_rows,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 1967 ierr = MatSetType(B,type);CHKERRQ(ierr); 1968 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(B,bs,0,s_browlengths);CHKERRQ(ierr); 1969 a = (Mat_SeqSBAIJ*)B->data; 1970 1971 /* set matrix "i" values */ 1972 a->i[0] = 0; 1973 for (i=1; i<= mbs; i++) { 1974 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 1975 a->ilen[i-1] = s_browlengths[i-1]; 1976 } 1977 a->nz = a->i[mbs]; 1978 1979 /* read in nonzero values */ 1980 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 1981 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 1982 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 1983 1984 /* set "a" and "j" values into matrix */ 1985 nzcount = 0; jcount = 0; 1986 for (i=0; i<mbs; i++) { 1987 nzcountb = nzcount; 1988 nmask = 0; 1989 for (j=0; j<bs; j++) { 1990 kmax = rowlengths[i*bs+j]; 1991 for (k=0; k<kmax; k++) { 1992 tmp = jj[nzcount++]/bs; /* block col. index */ 1993 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 1994 } 1995 } 1996 /* sort the masked values */ 1997 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 1998 1999 /* set "j" values into matrix */ 2000 maskcount = 1; 2001 for (j=0; j<nmask; j++) { 2002 a->j[jcount++] = masked[j]; 2003 mask[masked[j]] = maskcount++; 2004 } 2005 2006 /* set "a" values into matrix */ 2007 ishift = bs2*a->i[i]; 2008 for (j=0; j<bs; j++) { 2009 kmax = rowlengths[i*bs+j]; 2010 for (k=0; k<kmax; k++) { 2011 tmp = jj[nzcountb]/bs ; /* block col. index */ 2012 if (tmp >= i){ 2013 block = mask[tmp] - 1; 2014 point = jj[nzcountb] - bs*tmp; 2015 idx = ishift + bs2*block + j + bs*point; 2016 a->a[idx] = aa[nzcountb]; 2017 } 2018 nzcountb++; 2019 } 2020 } 2021 /* zero out the mask elements we set */ 2022 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2023 } 2024 if (jcount != a->nz) SETERRQ(PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2025 2026 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2027 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2028 ierr = PetscFree(aa);CHKERRQ(ierr); 2029 ierr = PetscFree(jj);CHKERRQ(ierr); 2030 ierr = PetscFree(mask);CHKERRQ(ierr); 2031 2032 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2033 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2034 ierr = MatView_Private(B);CHKERRQ(ierr); 2035 *A = B; 2036 PetscFunctionReturn(0); 2037 } 2038 2039 #undef __FUNCT__ 2040 #define __FUNCT__ "MatRelax_SeqSBAIJ" 2041 PetscErrorCode MatRelax_SeqSBAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 2042 { 2043 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 2044 MatScalar *aa=a->a,*v,*v1; 2045 PetscScalar *x,*b,*t,sum,d; 2046 PetscErrorCode ierr; 2047 PetscInt m=a->mbs,bs=A->rmap->bs,*ai=a->i,*aj=a->j; 2048 PetscInt nz,nz1,*vj,*vj1,i; 2049 2050 PetscFunctionBegin; 2051 if (flag & SOR_EISENSTAT) SETERRQ(PETSC_ERR_SUP,"No support yet for Eisenstat"); 2052 2053 its = its*lits; 2054 if (its <= 0) SETERRQ2(PETSC_ERR_ARG_WRONG,"Relaxation requires global its %D and local its %D both positive",its,lits); 2055 2056 if (bs > 1) 2057 SETERRQ(PETSC_ERR_SUP,"SSOR for block size > 1 is not yet implemented"); 2058 2059 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 2060 if (xx != bb) { 2061 ierr = VecGetArray(bb,&b);CHKERRQ(ierr); 2062 } else { 2063 b = x; 2064 } 2065 2066 if (!a->relax_work) { 2067 ierr = PetscMalloc(m*sizeof(PetscScalar),&a->relax_work);CHKERRQ(ierr); 2068 } 2069 t = a->relax_work; 2070 2071 if (flag & SOR_ZERO_INITIAL_GUESS) { 2072 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2073 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2074 2075 for (i=0; i<m; i++){ 2076 d = *(aa + ai[i]); /* diag[i] */ 2077 v = aa + ai[i] + 1; 2078 vj = aj + ai[i] + 1; 2079 nz = ai[i+1] - ai[i] - 1; 2080 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2081 x[i] = omega*t[i]/d; 2082 while (nz--) t[*vj++] -= x[i]*(*v++); /* update rhs */ 2083 } 2084 } 2085 2086 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2087 if (!(flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP)){ 2088 t = b; 2089 } 2090 2091 for (i=m-1; i>=0; i--){ 2092 d = *(aa + ai[i]); 2093 v = aa + ai[i] + 1; 2094 vj = aj + ai[i] + 1; 2095 nz = ai[i+1] - ai[i] - 1; 2096 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2097 sum = t[i]; 2098 while (nz--) sum -= x[*vj++]*(*v++); 2099 x[i] = (1-omega)*x[i] + omega*sum/d; 2100 } 2101 t = a->relax_work; 2102 } 2103 its--; 2104 } 2105 2106 while (its--) { 2107 /* 2108 forward sweep: 2109 for i=0,...,m-1: 2110 sum[i] = (b[i] - U(i,:)x )/d[i]; 2111 x[i] = (1-omega)x[i] + omega*sum[i]; 2112 b = b - x[i]*U^T(i,:); 2113 2114 */ 2115 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP){ 2116 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2117 2118 for (i=0; i<m; i++){ 2119 d = *(aa + ai[i]); /* diag[i] */ 2120 v = aa + ai[i] + 1; v1=v; 2121 vj = aj + ai[i] + 1; vj1=vj; 2122 nz = ai[i+1] - ai[i] - 1; nz1=nz; 2123 sum = t[i]; 2124 ierr = PetscLogFlops(4*nz-2);CHKERRQ(ierr); 2125 while (nz1--) sum -= (*v1++)*x[*vj1++]; 2126 x[i] = (1-omega)*x[i] + omega*sum/d; 2127 while (nz--) t[*vj++] -= x[i]*(*v++); 2128 } 2129 } 2130 2131 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP){ 2132 /* 2133 backward sweep: 2134 b = b - x[i]*U^T(i,:), i=0,...,n-2 2135 for i=m-1,...,0: 2136 sum[i] = (b[i] - U(i,:)x )/d[i]; 2137 x[i] = (1-omega)x[i] + omega*sum[i]; 2138 */ 2139 /* if there was a forward sweep done above then I thing the next two for loops are not needed */ 2140 ierr = PetscMemcpy(t,b,m*sizeof(PetscScalar));CHKERRQ(ierr); 2141 2142 for (i=0; i<m-1; i++){ /* update rhs */ 2143 v = aa + ai[i] + 1; 2144 vj = aj + ai[i] + 1; 2145 nz = ai[i+1] - ai[i] - 1; 2146 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2147 while (nz--) t[*vj++] -= x[i]*(*v++); 2148 } 2149 for (i=m-1; i>=0; i--){ 2150 d = *(aa + ai[i]); 2151 v = aa + ai[i] + 1; 2152 vj = aj + ai[i] + 1; 2153 nz = ai[i+1] - ai[i] - 1; 2154 ierr = PetscLogFlops(2*nz-1);CHKERRQ(ierr); 2155 sum = t[i]; 2156 while (nz--) sum -= x[*vj++]*(*v++); 2157 x[i] = (1-omega)*x[i] + omega*sum/d; 2158 } 2159 } 2160 } 2161 2162 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 2163 if (bb != xx) { 2164 ierr = VecRestoreArray(bb,&b);CHKERRQ(ierr); 2165 } 2166 PetscFunctionReturn(0); 2167 } 2168 2169 #undef __FUNCT__ 2170 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2171 /*@ 2172 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2173 (upper triangular entries in CSR format) provided by the user. 2174 2175 Collective on MPI_Comm 2176 2177 Input Parameters: 2178 + comm - must be an MPI communicator of size 1 2179 . bs - size of block 2180 . m - number of rows 2181 . n - number of columns 2182 . i - row indices 2183 . j - column indices 2184 - a - matrix values 2185 2186 Output Parameter: 2187 . mat - the matrix 2188 2189 Level: intermediate 2190 2191 Notes: 2192 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2193 once the matrix is destroyed 2194 2195 You cannot set new nonzero locations into this matrix, that will generate an error. 2196 2197 The i and j indices are 0 based 2198 2199 .seealso: MatCreate(), MatCreateMPISBAIJ(), MatCreateSeqSBAIJ() 2200 2201 @*/ 2202 PetscErrorCode PETSCMAT_DLLEXPORT MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 2203 { 2204 PetscErrorCode ierr; 2205 PetscInt ii; 2206 Mat_SeqSBAIJ *sbaij; 2207 2208 PetscFunctionBegin; 2209 if (bs != 1) SETERRQ1(PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2210 if (i[0]) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2211 2212 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2213 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2214 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2215 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2216 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2217 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2218 2219 sbaij->i = i; 2220 sbaij->j = j; 2221 sbaij->a = a; 2222 sbaij->singlemalloc = PETSC_FALSE; 2223 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2224 sbaij->free_a = PETSC_FALSE; 2225 sbaij->free_ij = PETSC_FALSE; 2226 2227 for (ii=0; ii<m; ii++) { 2228 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2229 #if defined(PETSC_USE_DEBUG) 2230 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]); 2231 #endif 2232 } 2233 #if defined(PETSC_USE_DEBUG) 2234 for (ii=0; ii<sbaij->i[m]; ii++) { 2235 if (j[ii] < 0) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2236 if (j[ii] > n - 1) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %d index = %d",ii,j[ii]); 2237 } 2238 #endif 2239 2240 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2241 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2242 PetscFunctionReturn(0); 2243 } 2244 2245 2246 2247 2248 2249