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