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