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