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