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