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