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