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