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