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