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((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,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 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 = PetscMalloc((*ncols+row)*sizeof(PetscScalar),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 = PetscMalloc((*ncols+row)*sizeof(PetscInt),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 static PetscErrorCode MatView_SeqSBAIJ_ASCII(Mat A,PetscViewer viewer) 350 { 351 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 352 PetscErrorCode ierr; 353 PetscInt i,j,bs = A->rmap->bs,k,l,bs2=a->bs2; 354 PetscViewerFormat format; 355 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 PetscRealPart(a->a[bs2*k + l*bs + j]),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 PetscRealPart(a->a[bs2*k + l*bs + j]),-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,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,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 ierr = PetscObjectPrintClassNamePrefixType((PetscObject)A,viewer);CHKERRQ(ierr); 403 if (A->factortype) { /* for factored matrix */ 404 if (bs>1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"matrix is factored with bs>1. Not implemented yet"); 405 406 diag=a->diag; 407 for (i=0; i<a->mbs; i++) { /* for row block i */ 408 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i);CHKERRQ(ierr); 409 /* diagonal entry */ 410 #if defined(PETSC_USE_COMPLEX) 411 if (PetscImaginaryPart(a->a[diag[i]]) > 0.0) { 412 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); 413 } else if (PetscImaginaryPart(a->a[diag[i]]) < 0.0) { 414 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); 415 } else { 416 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],PetscRealPart(1.0/a->a[diag[i]]));CHKERRQ(ierr); 417 } 418 #else 419 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[diag[i]],1.0/a->a[diag[i]]);CHKERRQ(ierr); 420 #endif 421 /* off-diagonal entries */ 422 for (k=a->i[i]; k<a->i[i+1]-1; k++) { 423 #if defined(PETSC_USE_COMPLEX) 424 if (PetscImaginaryPart(a->a[k]) > 0.0) { 425 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k],PetscRealPart(a->a[k]),PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 426 } else if (PetscImaginaryPart(a->a[k]) < 0.0) { 427 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k],PetscRealPart(a->a[k]),-PetscImaginaryPart(a->a[k]));CHKERRQ(ierr); 428 } else { 429 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k],PetscRealPart(a->a[k]));CHKERRQ(ierr); 430 } 431 #else 432 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",a->j[k],a->a[k]);CHKERRQ(ierr); 433 #endif 434 } 435 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 436 } 437 438 } else { /* for non-factored matrix */ 439 for (i=0; i<a->mbs; i++) { /* for row block i */ 440 for (j=0; j<bs; j++) { /* for row bs*i + j */ 441 ierr = PetscViewerASCIIPrintf(viewer,"row %D:",i*bs+j);CHKERRQ(ierr); 442 for (k=a->i[i]; k<a->i[i+1]; k++) { /* for column block */ 443 for (l=0; l<bs; l++) { /* for column */ 444 #if defined(PETSC_USE_COMPLEX) 445 if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) > 0.0) { 446 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G + %G i) ",bs*a->j[k]+l, 447 PetscRealPart(a->a[bs2*k + l*bs + j]),PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 448 } else if (PetscImaginaryPart(a->a[bs2*k + l*bs + j]) < 0.0) { 449 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G - %G i) ",bs*a->j[k]+l, 450 PetscRealPart(a->a[bs2*k + l*bs + j]),-PetscImaginaryPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 451 } else { 452 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,PetscRealPart(a->a[bs2*k + l*bs + j]));CHKERRQ(ierr); 453 } 454 #else 455 ierr = PetscViewerASCIIPrintf(viewer," (%D, %G) ",bs*a->j[k]+l,a->a[bs2*k + l*bs + j]);CHKERRQ(ierr); 456 #endif 457 } 458 } 459 ierr = PetscViewerASCIIPrintf(viewer,"\n");CHKERRQ(ierr); 460 } 461 } 462 } 463 ierr = PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);CHKERRQ(ierr); 464 } 465 ierr = PetscViewerFlush(viewer);CHKERRQ(ierr); 466 PetscFunctionReturn(0); 467 } 468 469 #include <petscdraw.h> 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",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) stepval = (n-1)*bs; 658 else stepval = (m-1)*bs; 659 660 for (k=0; k<m; k++) { /* loop over added rows */ 661 row = im[k]; 662 if (row < 0) continue; 663 #if defined(PETSC_USE_DEBUG) 664 if (row >= a->mbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %D max %D",row,a->mbs-1); 665 #endif 666 rp = aj + ai[row]; 667 ap = aa + bs2*ai[row]; 668 rmax = imax[row]; 669 nrow = ailen[row]; 670 low = 0; 671 high = nrow; 672 for (l=0; l<n; l++) { /* loop over added columns */ 673 if (in[l] < 0) continue; 674 col = in[l]; 675 #if defined(PETSC_USE_DEBUG) 676 if (col >= a->nbs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %D max %D",col,a->nbs-1); 677 #endif 678 if (col < row) { 679 if (a->ignore_ltriangular) continue; /* ignore lower triangular block */ 680 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)"); 681 } 682 if (roworiented) value = v + k*(stepval+bs)*bs + l*bs; 683 else value = v + l*(stepval+bs)*bs + k*bs; 684 685 if (col <= lastcol) low = 0; 686 else high = nrow; 687 688 lastcol = col; 689 while (high-low > 7) { 690 t = (low+high)/2; 691 if (rp[t] > col) high = t; 692 else low = t; 693 } 694 for (i=low; i<high; i++) { 695 if (rp[i] > col) break; 696 if (rp[i] == col) { 697 bap = ap + bs2*i; 698 if (roworiented) { 699 if (is == ADD_VALUES) { 700 for (ii=0; ii<bs; ii++,value+=stepval) { 701 for (jj=ii; jj<bs2; jj+=bs) { 702 bap[jj] += *value++; 703 } 704 } 705 } else { 706 for (ii=0; ii<bs; ii++,value+=stepval) { 707 for (jj=ii; jj<bs2; jj+=bs) { 708 bap[jj] = *value++; 709 } 710 } 711 } 712 } else { 713 if (is == ADD_VALUES) { 714 for (ii=0; ii<bs; ii++,value+=stepval) { 715 for (jj=0; jj<bs; jj++) { 716 *bap++ += *value++; 717 } 718 } 719 } else { 720 for (ii=0; ii<bs; ii++,value+=stepval) { 721 for (jj=0; jj<bs; jj++) { 722 *bap++ = *value++; 723 } 724 } 725 } 726 } 727 goto noinsert2; 728 } 729 } 730 if (nonew == 1) goto noinsert2; 731 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 732 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 733 N = nrow++ - 1; high++; 734 /* shift up all the later entries in this row */ 735 for (ii=N; ii>=i; ii--) { 736 rp[ii+1] = rp[ii]; 737 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 738 } 739 if (N >= i) { 740 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 741 } 742 rp[i] = col; 743 bap = ap + bs2*i; 744 if (roworiented) { 745 for (ii=0; ii<bs; ii++,value+=stepval) { 746 for (jj=ii; jj<bs2; jj+=bs) { 747 bap[jj] = *value++; 748 } 749 } 750 } else { 751 for (ii=0; ii<bs; ii++,value+=stepval) { 752 for (jj=0; jj<bs; jj++) { 753 *bap++ = *value++; 754 } 755 } 756 } 757 noinsert2:; 758 low = i; 759 } 760 ailen[row] = nrow; 761 } 762 PetscFunctionReturn(0); 763 } 764 765 /* 766 This is not yet used 767 */ 768 #undef __FUNCT__ 769 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode" 770 PetscErrorCode MatAssemblyEnd_SeqSBAIJ_SeqAIJ_Inode(Mat A) 771 { 772 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 773 PetscErrorCode ierr; 774 const PetscInt *ai = a->i, *aj = a->j,*cols; 775 PetscInt i = 0,j,blk_size,m = A->rmap->n,node_count = 0,nzx,nzy,*ns,row,nz,cnt,cnt2,*counts; 776 PetscBool flag; 777 778 PetscFunctionBegin; 779 ierr = PetscMalloc(m*sizeof(PetscInt),&ns);CHKERRQ(ierr); 780 while (i < m) { 781 nzx = ai[i+1] - ai[i]; /* Number of nonzeros */ 782 /* Limits the number of elements in a node to 'a->inode.limit' */ 783 for (j=i+1,blk_size=1; j<m && blk_size <a->inode.limit; ++j,++blk_size) { 784 nzy = ai[j+1] - ai[j]; 785 if (nzy != (nzx - j + i)) break; 786 ierr = PetscMemcmp(aj + ai[i] + j - i,aj + ai[j],nzy*sizeof(PetscInt),&flag);CHKERRQ(ierr); 787 if (!flag) break; 788 } 789 ns[node_count++] = blk_size; 790 791 i = j; 792 } 793 if (!a->inode.size && m && node_count > .9*m) { 794 ierr = PetscFree(ns);CHKERRQ(ierr); 795 ierr = PetscInfo2(A,"Found %D nodes out of %D rows. Not using Inode routines\n",node_count,m);CHKERRQ(ierr); 796 } else { 797 a->inode.node_count = node_count; 798 799 ierr = PetscMalloc(node_count*sizeof(PetscInt),&a->inode.size);CHKERRQ(ierr); 800 ierr = PetscLogObjectMemory((PetscObject)A,node_count*sizeof(PetscInt));CHKERRQ(ierr); 801 ierr = PetscMemcpy(a->inode.size,ns,node_count*sizeof(PetscInt));CHKERRQ(ierr); 802 ierr = PetscFree(ns);CHKERRQ(ierr); 803 ierr = PetscInfo3(A,"Found %D nodes of %D. Limit used: %D. Using Inode routines\n",node_count,m,a->inode.limit);CHKERRQ(ierr); 804 805 /* count collections of adjacent columns in each inode */ 806 row = 0; 807 cnt = 0; 808 for (i=0; i<node_count; i++) { 809 cols = aj + ai[row] + a->inode.size[i]; 810 nz = ai[row+1] - ai[row] - a->inode.size[i]; 811 for (j=1; j<nz; j++) { 812 if (cols[j] != cols[j-1]+1) cnt++; 813 } 814 cnt++; 815 row += a->inode.size[i]; 816 } 817 ierr = PetscMalloc(2*cnt*sizeof(PetscInt),&counts);CHKERRQ(ierr); 818 cnt = 0; 819 row = 0; 820 for (i=0; i<node_count; i++) { 821 cols = aj + ai[row] + a->inode.size[i]; 822 counts[2*cnt] = cols[0]; 823 nz = ai[row+1] - ai[row] - a->inode.size[i]; 824 cnt2 = 1; 825 for (j=1; j<nz; j++) { 826 if (cols[j] != cols[j-1]+1) { 827 counts[2*(cnt++)+1] = cnt2; 828 counts[2*cnt] = cols[j]; 829 cnt2 = 1; 830 } else cnt2++; 831 } 832 counts[2*(cnt++)+1] = cnt2; 833 row += a->inode.size[i]; 834 } 835 ierr = PetscIntView(2*cnt,counts,0);CHKERRQ(ierr); 836 } 837 PetscFunctionReturn(0); 838 } 839 840 #undef __FUNCT__ 841 #define __FUNCT__ "MatAssemblyEnd_SeqSBAIJ" 842 PetscErrorCode MatAssemblyEnd_SeqSBAIJ(Mat A,MatAssemblyType mode) 843 { 844 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 845 PetscErrorCode ierr; 846 PetscInt fshift = 0,i,j,*ai = a->i,*aj = a->j,*imax = a->imax; 847 PetscInt m = A->rmap->N,*ip,N,*ailen = a->ilen; 848 PetscInt mbs = a->mbs,bs2 = a->bs2,rmax = 0; 849 MatScalar *aa = a->a,*ap; 850 851 PetscFunctionBegin; 852 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 853 854 if (m) rmax = ailen[0]; 855 for (i=1; i<mbs; i++) { 856 /* move each row back by the amount of empty slots (fshift) before it*/ 857 fshift += imax[i-1] - ailen[i-1]; 858 rmax = PetscMax(rmax,ailen[i]); 859 if (fshift) { 860 ip = aj + ai[i]; ap = aa + bs2*ai[i]; 861 N = ailen[i]; 862 for (j=0; j<N; j++) { 863 ip[j-fshift] = ip[j]; 864 ierr = PetscMemcpy(ap+(j-fshift)*bs2,ap+j*bs2,bs2*sizeof(MatScalar));CHKERRQ(ierr); 865 } 866 } 867 ai[i] = ai[i-1] + ailen[i-1]; 868 } 869 if (mbs) { 870 fshift += imax[mbs-1] - ailen[mbs-1]; 871 ai[mbs] = ai[mbs-1] + ailen[mbs-1]; 872 } 873 /* reset ilen and imax for each row */ 874 for (i=0; i<mbs; i++) { 875 ailen[i] = imax[i] = ai[i+1] - ai[i]; 876 } 877 a->nz = ai[mbs]; 878 879 /* diagonals may have moved, reset it */ 880 if (a->diag) { 881 ierr = PetscMemcpy(a->diag,ai,mbs*sizeof(PetscInt));CHKERRQ(ierr); 882 } 883 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); 884 885 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); 886 ierr = PetscInfo1(A,"Number of mallocs during MatSetValues is %D\n",a->reallocs);CHKERRQ(ierr); 887 ierr = PetscInfo1(A,"Most nonzeros blocks in any row is %D\n",rmax);CHKERRQ(ierr); 888 889 A->info.mallocs += a->reallocs; 890 a->reallocs = 0; 891 A->info.nz_unneeded = (PetscReal)fshift*bs2; 892 a->idiagvalid = PETSC_FALSE; 893 894 if (A->cmap->n < 65536 && A->cmap->bs == 1) { 895 if (a->jshort && a->free_jshort) { 896 /* when matrix data structure is changed, previous jshort must be replaced */ 897 ierr = PetscFree(a->jshort);CHKERRQ(ierr); 898 } 899 ierr = PetscMalloc(a->i[A->rmap->n]*sizeof(unsigned short),&a->jshort);CHKERRQ(ierr); 900 ierr = PetscLogObjectMemory((PetscObject)A,a->i[A->rmap->n]*sizeof(unsigned short));CHKERRQ(ierr); 901 for (i=0; i<a->i[A->rmap->n]; i++) a->jshort[i] = a->j[i]; 902 A->ops->mult = MatMult_SeqSBAIJ_1_ushort; 903 A->ops->sor = MatSOR_SeqSBAIJ_ushort; 904 a->free_jshort = PETSC_TRUE; 905 } 906 PetscFunctionReturn(0); 907 } 908 909 /* 910 This function returns an array of flags which indicate the locations of contiguous 911 blocks that should be zeroed. for eg: if bs = 3 and is = [0,1,2,3,5,6,7,8,9] 912 then the resulting sizes = [3,1,1,3,1] correspondig to sets [(0,1,2),(3),(5),(6,7,8),(9)] 913 Assume: sizes should be long enough to hold all the values. 914 */ 915 #undef __FUNCT__ 916 #define __FUNCT__ "MatZeroRows_SeqSBAIJ_Check_Blocks" 917 PetscErrorCode MatZeroRows_SeqSBAIJ_Check_Blocks(PetscInt idx[],PetscInt n,PetscInt bs,PetscInt sizes[], PetscInt *bs_max) 918 { 919 PetscInt i,j,k,row; 920 PetscBool flg; 921 922 PetscFunctionBegin; 923 for (i=0,j=0; i<n; j++) { 924 row = idx[i]; 925 if (row%bs!=0) { /* Not the begining of a block */ 926 sizes[j] = 1; 927 i++; 928 } else if (i+bs > n) { /* Beginning of a block, but complete block doesn't exist (at idx end) */ 929 sizes[j] = 1; /* Also makes sure atleast 'bs' values exist for next else */ 930 i++; 931 } else { /* Begining of the block, so check if the complete block exists */ 932 flg = PETSC_TRUE; 933 for (k=1; k<bs; k++) { 934 if (row+k != idx[i+k]) { /* break in the block */ 935 flg = PETSC_FALSE; 936 break; 937 } 938 } 939 if (flg) { /* No break in the bs */ 940 sizes[j] = bs; 941 i += bs; 942 } else { 943 sizes[j] = 1; 944 i++; 945 } 946 } 947 } 948 *bs_max = j; 949 PetscFunctionReturn(0); 950 } 951 952 953 /* Only add/insert a(i,j) with i<=j (blocks). 954 Any a(i,j) with i>j input by user is ingored. 955 */ 956 957 #undef __FUNCT__ 958 #define __FUNCT__ "MatSetValues_SeqSBAIJ" 959 PetscErrorCode MatSetValues_SeqSBAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 960 { 961 Mat_SeqSBAIJ *a = (Mat_SeqSBAIJ*)A->data; 962 PetscErrorCode ierr; 963 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N,lastcol = -1; 964 PetscInt *imax=a->imax,*ai=a->i,*ailen=a->ilen,roworiented=a->roworiented; 965 PetscInt *aj =a->j,nonew=a->nonew,bs=A->rmap->bs,brow,bcol; 966 PetscInt ridx,cidx,bs2=a->bs2; 967 MatScalar *ap,value,*aa=a->a,*bap; 968 969 PetscFunctionBegin; 970 if (v) PetscValidScalarPointer(v,6); 971 for (k=0; k<m; k++) { /* loop over added rows */ 972 row = im[k]; /* row number */ 973 brow = row/bs; /* block row number */ 974 if (row < 0) continue; 975 #if defined(PETSC_USE_DEBUG) 976 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); 977 #endif 978 rp = aj + ai[brow]; /*ptr to beginning of column value of the row block*/ 979 ap = aa + bs2*ai[brow]; /*ptr to beginning of element value of the row block*/ 980 rmax = imax[brow]; /* maximum space allocated for this row */ 981 nrow = ailen[brow]; /* actual length of this row */ 982 low = 0; 983 984 for (l=0; l<n; l++) { /* loop over added columns */ 985 if (in[l] < 0) continue; 986 #if defined(PETSC_USE_DEBUG) 987 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); 988 #endif 989 col = in[l]; 990 bcol = col/bs; /* block col number */ 991 992 if (brow > bcol) { 993 if (a->ignore_ltriangular) continue; /* ignore lower triangular values */ 994 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)"); 995 } 996 997 ridx = row % bs; cidx = col % bs; /*row and col index inside the block */ 998 if ((brow==bcol && ridx<=cidx) || (brow<bcol)) { 999 /* element value a(k,l) */ 1000 if (roworiented) value = v[l + k*n]; 1001 else value = v[k + l*m]; 1002 1003 /* move pointer bap to a(k,l) quickly and add/insert value */ 1004 if (col <= lastcol) low = 0; 1005 high = nrow; 1006 lastcol = col; 1007 while (high-low > 7) { 1008 t = (low+high)/2; 1009 if (rp[t] > bcol) high = t; 1010 else low = t; 1011 } 1012 for (i=low; i<high; i++) { 1013 if (rp[i] > bcol) break; 1014 if (rp[i] == bcol) { 1015 bap = ap + bs2*i + bs*cidx + ridx; 1016 if (is == ADD_VALUES) *bap += value; 1017 else *bap = value; 1018 /* for diag block, add/insert its symmetric element a(cidx,ridx) */ 1019 if (brow == bcol && ridx < cidx) { 1020 bap = ap + bs2*i + bs*ridx + cidx; 1021 if (is == ADD_VALUES) *bap += value; 1022 else *bap = value; 1023 } 1024 goto noinsert1; 1025 } 1026 } 1027 1028 if (nonew == 1) goto noinsert1; 1029 if (nonew == -1) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero (%D, %D) in the matrix", row, col); 1030 MatSeqXAIJReallocateAIJ(A,a->mbs,bs2,nrow,brow,bcol,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 1031 1032 N = nrow++ - 1; high++; 1033 /* shift up all the later entries in this row */ 1034 for (ii=N; ii>=i; ii--) { 1035 rp[ii+1] = rp[ii]; 1036 ierr = PetscMemcpy(ap+bs2*(ii+1),ap+bs2*(ii),bs2*sizeof(MatScalar));CHKERRQ(ierr); 1037 } 1038 if (N>=i) { 1039 ierr = PetscMemzero(ap+bs2*i,bs2*sizeof(MatScalar));CHKERRQ(ierr); 1040 } 1041 rp[i] = bcol; 1042 ap[bs2*i + bs*cidx + ridx] = value; 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 = PetscMalloc((inA->rmap->N+inA->rmap->bs)*sizeof(PetscScalar),&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,(PetscReal)(bs2*x->nz)/(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 A->same_nonzero = PETSC_TRUE; 1315 1316 /* zero the columns */ 1317 ierr = PetscMalloc(A->rmap->n*sizeof(PetscBool),&zeroed);CHKERRQ(ierr); 1318 ierr = PetscMemzero(zeroed,A->rmap->n*sizeof(PetscBool));CHKERRQ(ierr); 1319 for (i=0; i<is_n; i++) { 1320 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]); 1321 zeroed[is_idx[i]] = PETSC_TRUE; 1322 } 1323 if (vecs) { 1324 for (i=0; i<A->rmap->N; i++) { 1325 row = i/bs; 1326 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1327 for (k=0; k<bs; k++) { 1328 col = bs*baij->j[j] + k; 1329 if (col <= i) continue; 1330 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1331 if (!zeroed[i] && zeroed[col]) bb[i] -= aa[0]*xx[col]; 1332 if (zeroed[i] && !zeroed[col]) bb[col] -= aa[0]*xx[i]; 1333 } 1334 } 1335 } 1336 for (i=0; i<is_n; i++) bb[is_idx[i]] = diag*xx[is_idx[i]]; 1337 } 1338 1339 for (i=0; i<A->rmap->N; i++) { 1340 if (!zeroed[i]) { 1341 row = i/bs; 1342 for (j=baij->i[row]; j<baij->i[row+1]; j++) { 1343 for (k=0; k<bs; k++) { 1344 col = bs*baij->j[j] + k; 1345 if (zeroed[col]) { 1346 aa = ((MatScalar*)(baij->a)) + j*bs2 + (i%bs) + bs*k; 1347 aa[0] = 0.0; 1348 } 1349 } 1350 } 1351 } 1352 } 1353 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1354 if (vecs) { 1355 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1356 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1357 } 1358 1359 /* zero the rows */ 1360 for (i=0; i<is_n; i++) { 1361 row = is_idx[i]; 1362 count = (baij->i[row/bs +1] - baij->i[row/bs])*bs; 1363 aa = ((MatScalar*)(baij->a)) + baij->i[row/bs]*bs2 + (row%bs); 1364 for (k=0; k<count; k++) { 1365 aa[0] = zero; 1366 aa += bs; 1367 } 1368 if (diag != 0.0) { 1369 ierr = (*A->ops->setvalues)(A,1,&row,1,&row,&diag,INSERT_VALUES);CHKERRQ(ierr); 1370 } 1371 } 1372 ierr = MatAssemblyEnd_SeqSBAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1373 PetscFunctionReturn(0); 1374 } 1375 1376 /* -------------------------------------------------------------------*/ 1377 static struct _MatOps MatOps_Values = {MatSetValues_SeqSBAIJ, 1378 MatGetRow_SeqSBAIJ, 1379 MatRestoreRow_SeqSBAIJ, 1380 MatMult_SeqSBAIJ_N, 1381 /* 4*/ MatMultAdd_SeqSBAIJ_N, 1382 MatMult_SeqSBAIJ_N, /* transpose versions are same as non-transpose versions */ 1383 MatMultAdd_SeqSBAIJ_N, 1384 0, 1385 0, 1386 0, 1387 /* 10*/ 0, 1388 0, 1389 MatCholeskyFactor_SeqSBAIJ, 1390 MatSOR_SeqSBAIJ, 1391 MatTranspose_SeqSBAIJ, 1392 /* 15*/ MatGetInfo_SeqSBAIJ, 1393 MatEqual_SeqSBAIJ, 1394 MatGetDiagonal_SeqSBAIJ, 1395 MatDiagonalScale_SeqSBAIJ, 1396 MatNorm_SeqSBAIJ, 1397 /* 20*/ 0, 1398 MatAssemblyEnd_SeqSBAIJ, 1399 MatSetOption_SeqSBAIJ, 1400 MatZeroEntries_SeqSBAIJ, 1401 /* 24*/ 0, 1402 0, 1403 0, 1404 0, 1405 0, 1406 /* 29*/ MatSetUp_SeqSBAIJ, 1407 0, 1408 0, 1409 0, 1410 0, 1411 /* 34*/ MatDuplicate_SeqSBAIJ, 1412 0, 1413 0, 1414 0, 1415 MatICCFactor_SeqSBAIJ, 1416 /* 39*/ MatAXPY_SeqSBAIJ, 1417 MatGetSubMatrices_SeqSBAIJ, 1418 MatIncreaseOverlap_SeqSBAIJ, 1419 MatGetValues_SeqSBAIJ, 1420 MatCopy_SeqSBAIJ, 1421 /* 44*/ 0, 1422 MatScale_SeqSBAIJ, 1423 0, 1424 0, 1425 MatZeroRowsColumns_SeqSBAIJ, 1426 /* 49*/ 0, 1427 MatGetRowIJ_SeqSBAIJ, 1428 MatRestoreRowIJ_SeqSBAIJ, 1429 0, 1430 0, 1431 /* 54*/ 0, 1432 0, 1433 0, 1434 0, 1435 MatSetValuesBlocked_SeqSBAIJ, 1436 /* 59*/ MatGetSubMatrix_SeqSBAIJ, 1437 0, 1438 0, 1439 0, 1440 0, 1441 /* 64*/ 0, 1442 0, 1443 0, 1444 0, 1445 0, 1446 /* 69*/ MatGetRowMaxAbs_SeqSBAIJ, 1447 0, 1448 0, 1449 0, 1450 0, 1451 /* 74*/ 0, 1452 0, 1453 0, 1454 0, 1455 0, 1456 /* 79*/ 0, 1457 0, 1458 0, 1459 MatGetInertia_SeqSBAIJ, 1460 MatLoad_SeqSBAIJ, 1461 /* 84*/ MatIsSymmetric_SeqSBAIJ, 1462 MatIsHermitian_SeqSBAIJ, 1463 MatIsStructurallySymmetric_SeqSBAIJ, 1464 0, 1465 0, 1466 /* 89*/ 0, 1467 0, 1468 0, 1469 0, 1470 0, 1471 /* 94*/ 0, 1472 0, 1473 0, 1474 0, 1475 0, 1476 /* 99*/ 0, 1477 0, 1478 0, 1479 0, 1480 0, 1481 /*104*/ 0, 1482 MatRealPart_SeqSBAIJ, 1483 MatImaginaryPart_SeqSBAIJ, 1484 MatGetRowUpperTriangular_SeqSBAIJ, 1485 MatRestoreRowUpperTriangular_SeqSBAIJ, 1486 /*109*/ 0, 1487 0, 1488 0, 1489 0, 1490 MatMissingDiagonal_SeqSBAIJ, 1491 /*114*/ 0, 1492 0, 1493 0, 1494 0, 1495 0, 1496 /*119*/ 0, 1497 0, 1498 0, 1499 0, 1500 0, 1501 /*124*/ 0, 1502 0, 1503 0, 1504 0, 1505 0, 1506 /*129*/ 0, 1507 0, 1508 0, 1509 0, 1510 0, 1511 /*134*/ 0, 1512 0, 1513 0, 1514 0, 1515 0, 1516 /*139*/ 0, 1517 0, 1518 0 1519 }; 1520 1521 #undef __FUNCT__ 1522 #define __FUNCT__ "MatStoreValues_SeqSBAIJ" 1523 PetscErrorCode MatStoreValues_SeqSBAIJ(Mat mat) 1524 { 1525 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1526 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1527 PetscErrorCode ierr; 1528 1529 PetscFunctionBegin; 1530 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1531 1532 /* allocate space for values if not already there */ 1533 if (!aij->saved_values) { 1534 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 1535 } 1536 1537 /* copy values over */ 1538 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1539 PetscFunctionReturn(0); 1540 } 1541 1542 #undef __FUNCT__ 1543 #define __FUNCT__ "MatRetrieveValues_SeqSBAIJ" 1544 PetscErrorCode MatRetrieveValues_SeqSBAIJ(Mat mat) 1545 { 1546 Mat_SeqSBAIJ *aij = (Mat_SeqSBAIJ*)mat->data; 1547 PetscErrorCode ierr; 1548 PetscInt nz = aij->i[mat->rmap->N]*mat->rmap->bs*aij->bs2; 1549 1550 PetscFunctionBegin; 1551 if (aij->nonew != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 1552 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 1553 1554 /* copy values over */ 1555 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 1556 PetscFunctionReturn(0); 1557 } 1558 1559 #undef __FUNCT__ 1560 #define __FUNCT__ "MatSeqSBAIJSetPreallocation_SeqSBAIJ" 1561 PetscErrorCode MatSeqSBAIJSetPreallocation_SeqSBAIJ(Mat B,PetscInt bs,PetscInt nz,PetscInt *nnz) 1562 { 1563 Mat_SeqSBAIJ *b = (Mat_SeqSBAIJ*)B->data; 1564 PetscErrorCode ierr; 1565 PetscInt i,mbs,bs2; 1566 PetscBool skipallocation = PETSC_FALSE,flg = PETSC_FALSE,realalloc = PETSC_FALSE; 1567 1568 PetscFunctionBegin; 1569 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 1570 B->preallocated = PETSC_TRUE; 1571 1572 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1573 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1574 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1575 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1576 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1577 1578 mbs = B->rmap->N/bs; 1579 bs2 = bs*bs; 1580 1581 if (mbs*bs != B->rmap->N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Number rows, cols must be divisible by blocksize"); 1582 1583 if (nz == MAT_SKIP_ALLOCATION) { 1584 skipallocation = PETSC_TRUE; 1585 nz = 0; 1586 } 1587 1588 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 3; 1589 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 1590 if (nnz) { 1591 for (i=0; i<mbs; i++) { 1592 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]); 1593 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); 1594 } 1595 } 1596 1597 B->ops->mult = MatMult_SeqSBAIJ_N; 1598 B->ops->multadd = MatMultAdd_SeqSBAIJ_N; 1599 B->ops->multtranspose = MatMult_SeqSBAIJ_N; 1600 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_N; 1601 1602 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1603 if (!flg) { 1604 switch (bs) { 1605 case 1: 1606 B->ops->mult = MatMult_SeqSBAIJ_1; 1607 B->ops->multadd = MatMultAdd_SeqSBAIJ_1; 1608 B->ops->multtranspose = MatMult_SeqSBAIJ_1; 1609 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_1; 1610 break; 1611 case 2: 1612 B->ops->mult = MatMult_SeqSBAIJ_2; 1613 B->ops->multadd = MatMultAdd_SeqSBAIJ_2; 1614 B->ops->multtranspose = MatMult_SeqSBAIJ_2; 1615 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_2; 1616 break; 1617 case 3: 1618 B->ops->mult = MatMult_SeqSBAIJ_3; 1619 B->ops->multadd = MatMultAdd_SeqSBAIJ_3; 1620 B->ops->multtranspose = MatMult_SeqSBAIJ_3; 1621 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_3; 1622 break; 1623 case 4: 1624 B->ops->mult = MatMult_SeqSBAIJ_4; 1625 B->ops->multadd = MatMultAdd_SeqSBAIJ_4; 1626 B->ops->multtranspose = MatMult_SeqSBAIJ_4; 1627 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_4; 1628 break; 1629 case 5: 1630 B->ops->mult = MatMult_SeqSBAIJ_5; 1631 B->ops->multadd = MatMultAdd_SeqSBAIJ_5; 1632 B->ops->multtranspose = MatMult_SeqSBAIJ_5; 1633 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_5; 1634 break; 1635 case 6: 1636 B->ops->mult = MatMult_SeqSBAIJ_6; 1637 B->ops->multadd = MatMultAdd_SeqSBAIJ_6; 1638 B->ops->multtranspose = MatMult_SeqSBAIJ_6; 1639 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_6; 1640 break; 1641 case 7: 1642 B->ops->mult = MatMult_SeqSBAIJ_7; 1643 B->ops->multadd = MatMultAdd_SeqSBAIJ_7; 1644 B->ops->multtranspose = MatMult_SeqSBAIJ_7; 1645 B->ops->multtransposeadd = MatMultAdd_SeqSBAIJ_7; 1646 break; 1647 } 1648 } 1649 1650 b->mbs = mbs; 1651 b->nbs = mbs; 1652 if (!skipallocation) { 1653 if (!b->imax) { 1654 ierr = PetscMalloc2(mbs,PetscInt,&b->imax,mbs,PetscInt,&b->ilen);CHKERRQ(ierr); 1655 1656 b->free_imax_ilen = PETSC_TRUE; 1657 1658 ierr = PetscLogObjectMemory((PetscObject)B,2*mbs*sizeof(PetscInt));CHKERRQ(ierr); 1659 } 1660 if (!nnz) { 1661 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 1662 else if (nz <= 0) nz = 1; 1663 for (i=0; i<mbs; i++) b->imax[i] = nz; 1664 nz = nz*mbs; /* total nz */ 1665 } else { 1666 nz = 0; 1667 for (i=0; i<mbs; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 1668 } 1669 /* b->ilen will count nonzeros in each block row so far. */ 1670 for (i=0; i<mbs; i++) b->ilen[i] = 0; 1671 /* nz=(nz+mbs)/2; */ /* total diagonal and superdiagonal nonzero blocks */ 1672 1673 /* allocate the matrix space */ 1674 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 1675 ierr = PetscMalloc3(bs2*nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->N+1,PetscInt,&b->i);CHKERRQ(ierr); 1676 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->N+1)*sizeof(PetscInt)+nz*(bs2*sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 1677 ierr = PetscMemzero(b->a,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 1678 ierr = PetscMemzero(b->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 1679 1680 b->singlemalloc = PETSC_TRUE; 1681 1682 /* pointer to beginning of each row */ 1683 b->i[0] = 0; 1684 for (i=1; i<mbs+1; i++) b->i[i] = b->i[i-1] + b->imax[i-1]; 1685 1686 b->free_a = PETSC_TRUE; 1687 b->free_ij = PETSC_TRUE; 1688 } else { 1689 b->free_a = PETSC_FALSE; 1690 b->free_ij = PETSC_FALSE; 1691 } 1692 1693 B->rmap->bs = bs; 1694 b->bs2 = bs2; 1695 b->nz = 0; 1696 b->maxnz = nz; 1697 1698 b->inew = 0; 1699 b->jnew = 0; 1700 b->anew = 0; 1701 b->a2anew = 0; 1702 b->permute = PETSC_FALSE; 1703 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 1704 PetscFunctionReturn(0); 1705 } 1706 1707 #undef __FUNCT__ 1708 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ" 1709 PetscErrorCode MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ(Mat B,PetscInt bs,const PetscInt ii[],const PetscInt jj[], const PetscScalar V[]) 1710 { 1711 PetscInt i,j,m,nz,nz_max=0,*nnz; 1712 PetscScalar *values=0; 1713 PetscBool roworiented = ((Mat_SeqSBAIJ*)B->data)->roworiented; 1714 PetscErrorCode ierr; 1715 PetscFunctionBegin; 1716 if (bs < 1) SETERRQ1(PetscObjectComm((PetscObject)B),PETSC_ERR_ARG_OUTOFRANGE,"Invalid block size specified, must be positive but it is %D",bs); 1717 ierr = PetscLayoutSetBlockSize(B->rmap,bs);CHKERRQ(ierr); 1718 ierr = PetscLayoutSetBlockSize(B->cmap,bs);CHKERRQ(ierr); 1719 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 1720 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 1721 ierr = PetscLayoutGetBlockSize(B->rmap,&bs);CHKERRQ(ierr); 1722 m = B->rmap->n/bs; 1723 1724 if (ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"ii[0] must be 0 but it is %D",ii[0]); 1725 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 1726 for (i=0; i<m; i++) { 1727 nz = ii[i+1] - ii[i]; 1728 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D has a negative number of columns %D",i,nz); 1729 nz_max = PetscMax(nz_max,nz); 1730 nnz[i] = nz; 1731 } 1732 ierr = MatSeqSBAIJSetPreallocation(B,bs,0,nnz);CHKERRQ(ierr); 1733 ierr = PetscFree(nnz);CHKERRQ(ierr); 1734 1735 values = (PetscScalar*)V; 1736 if (!values) { 1737 ierr = PetscMalloc(bs*bs*nz_max*sizeof(PetscScalar),&values);CHKERRQ(ierr); 1738 ierr = PetscMemzero(values,bs*bs*nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 1739 } 1740 for (i=0; i<m; i++) { 1741 PetscInt ncols = ii[i+1] - ii[i]; 1742 const PetscInt *icols = jj + ii[i]; 1743 if (!roworiented || bs == 1) { 1744 const PetscScalar *svals = values + (V ? (bs*bs*ii[i]) : 0); 1745 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,ncols,icols,svals,INSERT_VALUES);CHKERRQ(ierr); 1746 } else { 1747 for (j=0; j<ncols; j++) { 1748 const PetscScalar *svals = values + (V ? (bs*bs*(ii[i]+j)) : 0); 1749 ierr = MatSetValuesBlocked_SeqSBAIJ(B,1,&i,1,&icols[j],svals,INSERT_VALUES);CHKERRQ(ierr); 1750 } 1751 } 1752 } 1753 if (!V) { ierr = PetscFree(values);CHKERRQ(ierr); } 1754 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1755 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1756 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 1757 PetscFunctionReturn(0); 1758 } 1759 1760 /* 1761 This is used to set the numeric factorization for both Cholesky and ICC symbolic factorization 1762 */ 1763 #undef __FUNCT__ 1764 #define __FUNCT__ "MatSeqSBAIJSetNumericFactorization_inplace" 1765 PetscErrorCode MatSeqSBAIJSetNumericFactorization_inplace(Mat B,PetscBool natural) 1766 { 1767 PetscErrorCode ierr; 1768 PetscBool flg = PETSC_FALSE; 1769 PetscInt bs = B->rmap->bs; 1770 1771 PetscFunctionBegin; 1772 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_no_unroll",&flg,NULL);CHKERRQ(ierr); 1773 if (flg) bs = 8; 1774 1775 if (!natural) { 1776 switch (bs) { 1777 case 1: 1778 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_inplace; 1779 break; 1780 case 2: 1781 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2; 1782 break; 1783 case 3: 1784 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3; 1785 break; 1786 case 4: 1787 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4; 1788 break; 1789 case 5: 1790 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5; 1791 break; 1792 case 6: 1793 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6; 1794 break; 1795 case 7: 1796 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7; 1797 break; 1798 default: 1799 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N; 1800 break; 1801 } 1802 } else { 1803 switch (bs) { 1804 case 1: 1805 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_1_NaturalOrdering_inplace; 1806 break; 1807 case 2: 1808 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_2_NaturalOrdering; 1809 break; 1810 case 3: 1811 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_3_NaturalOrdering; 1812 break; 1813 case 4: 1814 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_4_NaturalOrdering; 1815 break; 1816 case 5: 1817 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_5_NaturalOrdering; 1818 break; 1819 case 6: 1820 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_6_NaturalOrdering; 1821 break; 1822 case 7: 1823 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_7_NaturalOrdering; 1824 break; 1825 default: 1826 B->ops->choleskyfactornumeric = MatCholeskyFactorNumeric_SeqSBAIJ_N_NaturalOrdering; 1827 break; 1828 } 1829 } 1830 PetscFunctionReturn(0); 1831 } 1832 1833 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqAIJ(Mat, MatType,MatReuse,Mat*); 1834 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqBAIJ(Mat, MatType,MatReuse,Mat*); 1835 1836 #undef __FUNCT__ 1837 #define __FUNCT__ "MatGetFactor_seqsbaij_petsc" 1838 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_petsc(Mat A,MatFactorType ftype,Mat *B) 1839 { 1840 PetscInt n = A->rmap->n; 1841 PetscErrorCode ierr; 1842 1843 PetscFunctionBegin; 1844 #if defined(PETSC_USE_COMPLEX) 1845 if (A->hermitian) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Hermitian Factor is not supported"); 1846 #endif 1847 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 1848 ierr = MatSetSizes(*B,n,n,n,n);CHKERRQ(ierr); 1849 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1850 ierr = MatSetType(*B,MATSEQSBAIJ);CHKERRQ(ierr); 1851 ierr = MatSeqSBAIJSetPreallocation(*B,A->rmap->bs,MAT_SKIP_ALLOCATION,NULL);CHKERRQ(ierr); 1852 1853 (*B)->ops->choleskyfactorsymbolic = MatCholeskyFactorSymbolic_SeqSBAIJ; 1854 (*B)->ops->iccfactorsymbolic = MatICCFactorSymbolic_SeqSBAIJ; 1855 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Factor type not supported"); 1856 (*B)->factortype = ftype; 1857 PetscFunctionReturn(0); 1858 } 1859 1860 #undef __FUNCT__ 1861 #define __FUNCT__ "MatGetFactorAvailable_seqsbaij_petsc" 1862 PetscErrorCode MatGetFactorAvailable_seqsbaij_petsc(Mat A,MatFactorType ftype,PetscBool *flg) 1863 { 1864 PetscFunctionBegin; 1865 if (ftype == MAT_FACTOR_CHOLESKY || ftype == MAT_FACTOR_ICC) { 1866 *flg = PETSC_TRUE; 1867 } else { 1868 *flg = PETSC_FALSE; 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 #if defined(PETSC_HAVE_MUMPS) 1874 PETSC_EXTERN PetscErrorCode MatGetFactor_sbaij_mumps(Mat,MatFactorType,Mat*); 1875 #endif 1876 #if defined(PETSC_HAVE_PASTIX) 1877 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_pastix(Mat,MatFactorType,Mat*); 1878 #endif 1879 #if defined(PETSC_HAVE_CHOLMOD) 1880 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_cholmod(Mat,MatFactorType,Mat*); 1881 #endif 1882 PETSC_EXTERN PetscErrorCode MatGetFactor_seqsbaij_sbstrm(Mat,MatFactorType,Mat*); 1883 1884 /*MC 1885 MATSEQSBAIJ - MATSEQSBAIJ = "seqsbaij" - A matrix type to be used for sequential symmetric block sparse matrices, 1886 based on block compressed sparse row format. Only the upper triangular portion of the matrix is stored. 1887 1888 For complex numbers by default this matrix is symmetric, NOT Hermitian symmetric. To make it Hermitian symmetric you 1889 can call MatSetOption(Mat, MAT_HERMITIAN); after MatAssemblyEnd() 1890 1891 Options Database Keys: 1892 . -mat_type seqsbaij - sets the matrix type to "seqsbaij" during a call to MatSetFromOptions() 1893 1894 Notes: By default if you insert values into the lower triangular part of the matrix they are simply ignored (since they are not 1895 stored and it is assumed they symmetric to the upper triangular). If you call MatSetOption(Mat,MAT_IGNORE_LOWER_TRIANGULAR,PETSC_FALSE) or use 1896 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. 1897 1898 1899 Level: beginner 1900 1901 .seealso: MatCreateSeqSBAIJ 1902 M*/ 1903 1904 PETSC_EXTERN PetscErrorCode MatConvert_SeqSBAIJ_SeqSBSTRM(Mat, MatType,MatReuse,Mat*); 1905 1906 #undef __FUNCT__ 1907 #define __FUNCT__ "MatCreate_SeqSBAIJ" 1908 PETSC_EXTERN PetscErrorCode MatCreate_SeqSBAIJ(Mat B) 1909 { 1910 Mat_SeqSBAIJ *b; 1911 PetscErrorCode ierr; 1912 PetscMPIInt size; 1913 PetscBool no_unroll = PETSC_FALSE,no_inode = PETSC_FALSE; 1914 1915 PetscFunctionBegin; 1916 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 1917 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Comm must be of size 1"); 1918 1919 ierr = PetscNewLog(B,Mat_SeqSBAIJ,&b);CHKERRQ(ierr); 1920 B->data = (void*)b; 1921 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 1922 1923 B->ops->destroy = MatDestroy_SeqSBAIJ; 1924 B->ops->view = MatView_SeqSBAIJ; 1925 b->row = 0; 1926 b->icol = 0; 1927 b->reallocs = 0; 1928 b->saved_values = 0; 1929 b->inode.limit = 5; 1930 b->inode.max_limit = 5; 1931 1932 b->roworiented = PETSC_TRUE; 1933 b->nonew = 0; 1934 b->diag = 0; 1935 b->solve_work = 0; 1936 b->mult_work = 0; 1937 B->spptr = 0; 1938 B->info.nz_unneeded = (PetscReal)b->maxnz*b->bs2; 1939 b->keepnonzeropattern = PETSC_FALSE; 1940 b->xtoy = 0; 1941 b->XtoY = 0; 1942 1943 b->inew = 0; 1944 b->jnew = 0; 1945 b->anew = 0; 1946 b->a2anew = 0; 1947 b->permute = PETSC_FALSE; 1948 1949 b->ignore_ltriangular = PETSC_TRUE; 1950 1951 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_ignore_lower_triangular",&b->ignore_ltriangular,NULL);CHKERRQ(ierr); 1952 1953 b->getrow_utriangular = PETSC_FALSE; 1954 1955 ierr = PetscOptionsGetBool(((PetscObject)B)->prefix,"-mat_getrow_uppertriangular",&b->getrow_utriangular,NULL);CHKERRQ(ierr); 1956 1957 #if defined(PETSC_HAVE_PASTIX) 1958 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_pastix_C",MatGetFactor_seqsbaij_pastix);CHKERRQ(ierr); 1959 #endif 1960 #if defined(PETSC_HAVE_MUMPS) 1961 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_mumps_C",MatGetFactor_sbaij_mumps);CHKERRQ(ierr); 1962 #endif 1963 #if defined(PETSC_HAVE_CHOLMOD) 1964 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_cholmod_C",MatGetFactor_seqsbaij_cholmod);CHKERRQ(ierr); 1965 #endif 1966 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactorAvailable_petsc_C",MatGetFactorAvailable_seqsbaij_petsc);CHKERRQ(ierr); 1967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_petsc_C",MatGetFactor_seqsbaij_petsc);CHKERRQ(ierr); 1968 ierr = PetscObjectComposeFunction((PetscObject)B,"MatGetFactor_sbstrm_C",MatGetFactor_seqsbaij_sbstrm);CHKERRQ(ierr); 1969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqSBAIJ);CHKERRQ(ierr); 1970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqSBAIJ);CHKERRQ(ierr); 1971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetColumnIndices_C",MatSeqSBAIJSetColumnIndices_SeqSBAIJ);CHKERRQ(ierr); 1972 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqaij_C",MatConvert_SeqSBAIJ_SeqAIJ);CHKERRQ(ierr); 1973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqbaij_C",MatConvert_SeqSBAIJ_SeqBAIJ);CHKERRQ(ierr); 1974 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocation_C",MatSeqSBAIJSetPreallocation_SeqSBAIJ);CHKERRQ(ierr); 1975 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqSBAIJSetPreallocationCSR_C",MatSeqSBAIJSetPreallocationCSR_SeqSBAIJ);CHKERRQ(ierr); 1976 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqsbaij_seqsbstrm_C",MatConvert_SeqSBAIJ_SeqSBSTRM);CHKERRQ(ierr); 1977 1978 B->symmetric = PETSC_TRUE; 1979 B->structurally_symmetric = PETSC_TRUE; 1980 B->symmetric_set = PETSC_TRUE; 1981 B->structurally_symmetric_set = PETSC_TRUE; 1982 1983 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQSBAIJ);CHKERRQ(ierr); 1984 1985 ierr = PetscOptionsBegin(PetscObjectComm((PetscObject)B),((PetscObject)B)->prefix,"Options for SEQSBAIJ matrix","Mat");CHKERRQ(ierr); 1986 ierr = PetscOptionsBool("-mat_no_unroll","Do not optimize for inodes (slower)",NULL,no_unroll,&no_unroll,NULL);CHKERRQ(ierr); 1987 if (no_unroll) { 1988 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_unroll\n");CHKERRQ(ierr); 1989 } 1990 ierr = PetscOptionsBool("-mat_no_inode","Do not optimize for inodes (slower)",NULL,no_inode,&no_inode,NULL);CHKERRQ(ierr); 1991 if (no_inode) { 1992 ierr = PetscInfo(B,"Not using Inode routines due to -mat_no_inode\n");CHKERRQ(ierr); 1993 } 1994 ierr = PetscOptionsInt("-mat_inode_limit","Do not use inodes larger then this value",NULL,b->inode.limit,&b->inode.limit,NULL);CHKERRQ(ierr); 1995 ierr = PetscOptionsEnd();CHKERRQ(ierr); 1996 b->inode.use = (PetscBool)(!(no_unroll || no_inode)); 1997 if (b->inode.limit > b->inode.max_limit) b->inode.limit = b->inode.max_limit; 1998 PetscFunctionReturn(0); 1999 } 2000 2001 #undef __FUNCT__ 2002 #define __FUNCT__ "MatSeqSBAIJSetPreallocation" 2003 /*@C 2004 MatSeqSBAIJSetPreallocation - Creates a sparse symmetric matrix in block AIJ (block 2005 compressed row) format. For good matrix assembly performance the 2006 user should preallocate the matrix storage by setting the parameter nz 2007 (or the array nnz). By setting these parameters accurately, performance 2008 during matrix assembly can be increased by more than a factor of 50. 2009 2010 Collective on Mat 2011 2012 Input Parameters: 2013 + A - the symmetric matrix 2014 . bs - size of block 2015 . nz - number of block nonzeros per block row (same for all rows) 2016 - nnz - array containing the number of block nonzeros in the upper triangular plus 2017 diagonal portion of each block (possibly different for each block row) or NULL 2018 2019 Options Database Keys: 2020 . -mat_no_unroll - uses code that does not unroll the loops in the 2021 block calculations (much slower) 2022 . -mat_block_size - size of the blocks to use (only works if a negative bs is passed in 2023 2024 Level: intermediate 2025 2026 Notes: 2027 Specify the preallocated storage with either nz or nnz (not both). 2028 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2029 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2030 2031 You can call MatGetInfo() to get information on how effective the preallocation was; 2032 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 2033 You can also run with the option -info and look for messages with the string 2034 malloc in them to see if additional memory allocation was needed. 2035 2036 If the nnz parameter is given then the nz parameter is ignored 2037 2038 2039 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2040 @*/ 2041 PetscErrorCode MatSeqSBAIJSetPreallocation(Mat B,PetscInt bs,PetscInt nz,const PetscInt nnz[]) 2042 { 2043 PetscErrorCode ierr; 2044 2045 PetscFunctionBegin; 2046 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2047 PetscValidType(B,1); 2048 PetscValidLogicalCollectiveInt(B,bs,2); 2049 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocation_C",(Mat,PetscInt,PetscInt,const PetscInt[]),(B,bs,nz,nnz));CHKERRQ(ierr); 2050 PetscFunctionReturn(0); 2051 } 2052 2053 #undef __FUNCT__ 2054 #define __FUNCT__ "MatSeqSBAIJSetPreallocationCSR" 2055 /*@C 2056 MatSeqSBAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in symmetric block AIJ format. 2057 2058 Input Parameters: 2059 + A - the matrix 2060 . i - the indices into j for the start of each local row (starts with zero) 2061 . j - the column indices for each local row (starts with zero) these must be sorted for each row 2062 - v - optional values in the matrix 2063 2064 Level: developer 2065 2066 Notes: 2067 The order of the entries in values is specified by the MatOption MAT_ROW_ORIENTED. For example, C programs 2068 may want to use the default MAT_ROW_ORIENTED=PETSC_TRUE and use an array v[nnz][bs][bs] where the second index is 2069 over rows within a block and the last index is over columns within a block row. Fortran programs will likely set 2070 MAT_ROW_ORIENTED=PETSC_FALSE and use a Fortran array v(bs,bs,nnz) in which the first index is over rows within a 2071 block column and the second index is over columns within a block. 2072 2073 .keywords: matrix, block, aij, compressed row, sparse 2074 2075 .seealso: MatCreate(), MatCreateSeqSBAIJ(), MatSetValuesBlocked(), MatSeqSBAIJSetPreallocation(), MATSEQSBAIJ 2076 @*/ 2077 PetscErrorCode MatSeqSBAIJSetPreallocationCSR(Mat B,PetscInt bs,const PetscInt i[],const PetscInt j[], const PetscScalar v[]) 2078 { 2079 PetscErrorCode ierr; 2080 2081 PetscFunctionBegin; 2082 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 2083 PetscValidType(B,1); 2084 PetscValidLogicalCollectiveInt(B,bs,2); 2085 ierr = PetscTryMethod(B,"MatSeqSBAIJSetPreallocationCSR_C",(Mat,PetscInt,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,bs,i,j,v));CHKERRQ(ierr); 2086 PetscFunctionReturn(0); 2087 } 2088 2089 #undef __FUNCT__ 2090 #define __FUNCT__ "MatCreateSeqSBAIJ" 2091 /*@C 2092 MatCreateSeqSBAIJ - Creates a sparse symmetric matrix in block AIJ (block 2093 compressed row) format. For good matrix assembly performance the 2094 user should preallocate the matrix storage by setting the parameter nz 2095 (or the array nnz). By setting these parameters accurately, performance 2096 during matrix assembly can be increased by more than a factor of 50. 2097 2098 Collective on MPI_Comm 2099 2100 Input Parameters: 2101 + comm - MPI communicator, set to PETSC_COMM_SELF 2102 . bs - size of block 2103 . m - number of rows, or number of columns 2104 . nz - number of block nonzeros per block row (same for all rows) 2105 - nnz - array containing the number of block nonzeros in the upper triangular plus 2106 diagonal portion of each block (possibly different for each block row) or NULL 2107 2108 Output Parameter: 2109 . A - the symmetric matrix 2110 2111 Options Database Keys: 2112 . -mat_no_unroll - uses code that does not unroll the loops in the 2113 block calculations (much slower) 2114 . -mat_block_size - size of the blocks to use 2115 2116 Level: intermediate 2117 2118 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 2119 MatXXXXSetPreallocation() paradgm instead of this routine directly. 2120 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 2121 2122 Notes: 2123 The number of rows and columns must be divisible by blocksize. 2124 This matrix type does not support complex Hermitian operation. 2125 2126 Specify the preallocated storage with either nz or nnz (not both). 2127 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 2128 allocation. See the <a href="../../docs/manual.pdf#nameddest=ch_mat">Mat chapter of the users manual</a> for details. 2129 2130 If the nnz parameter is given then the nz parameter is ignored 2131 2132 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatCreateSBAIJ() 2133 @*/ 2134 PetscErrorCode MatCreateSeqSBAIJ(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 2135 { 2136 PetscErrorCode ierr; 2137 2138 PetscFunctionBegin; 2139 ierr = MatCreate(comm,A);CHKERRQ(ierr); 2140 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 2141 ierr = MatSetType(*A,MATSEQSBAIJ);CHKERRQ(ierr); 2142 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*A,bs,nz,(PetscInt*)nnz);CHKERRQ(ierr); 2143 PetscFunctionReturn(0); 2144 } 2145 2146 #undef __FUNCT__ 2147 #define __FUNCT__ "MatDuplicate_SeqSBAIJ" 2148 PetscErrorCode MatDuplicate_SeqSBAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 2149 { 2150 Mat C; 2151 Mat_SeqSBAIJ *c,*a = (Mat_SeqSBAIJ*)A->data; 2152 PetscErrorCode ierr; 2153 PetscInt i,mbs = a->mbs,nz = a->nz,bs2 =a->bs2; 2154 2155 PetscFunctionBegin; 2156 if (a->i[mbs] != nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Corrupt matrix"); 2157 2158 *B = 0; 2159 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2160 ierr = MatSetSizes(C,A->rmap->N,A->cmap->n,A->rmap->N,A->cmap->n);CHKERRQ(ierr); 2161 ierr = MatSetType(C,MATSEQSBAIJ);CHKERRQ(ierr); 2162 ierr = PetscMemcpy(C->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr); 2163 c = (Mat_SeqSBAIJ*)C->data; 2164 2165 C->preallocated = PETSC_TRUE; 2166 C->factortype = A->factortype; 2167 c->row = 0; 2168 c->icol = 0; 2169 c->saved_values = 0; 2170 c->keepnonzeropattern = a->keepnonzeropattern; 2171 C->assembled = PETSC_TRUE; 2172 2173 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 2174 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 2175 c->bs2 = a->bs2; 2176 c->mbs = a->mbs; 2177 c->nbs = a->nbs; 2178 2179 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2180 c->imax = a->imax; 2181 c->ilen = a->ilen; 2182 c->free_imax_ilen = PETSC_FALSE; 2183 } else { 2184 ierr = PetscMalloc2((mbs+1),PetscInt,&c->imax,(mbs+1),PetscInt,&c->ilen);CHKERRQ(ierr); 2185 ierr = PetscLogObjectMemory((PetscObject)C,2*(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2186 for (i=0; i<mbs; i++) { 2187 c->imax[i] = a->imax[i]; 2188 c->ilen[i] = a->ilen[i]; 2189 } 2190 c->free_imax_ilen = PETSC_TRUE; 2191 } 2192 2193 /* allocate the matrix space */ 2194 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2195 ierr = PetscMalloc(bs2*nz*sizeof(MatScalar),&c->a);CHKERRQ(ierr); 2196 ierr = PetscLogObjectMemory((PetscObject)C,nz*bs2*sizeof(MatScalar));CHKERRQ(ierr); 2197 c->i = a->i; 2198 c->j = a->j; 2199 c->singlemalloc = PETSC_FALSE; 2200 c->free_a = PETSC_TRUE; 2201 c->free_ij = PETSC_FALSE; 2202 c->parent = A; 2203 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr); 2204 ierr = MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2205 ierr = MatSetOption(C,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 2206 } else { 2207 ierr = PetscMalloc3(bs2*nz,MatScalar,&c->a,nz,PetscInt,&c->j,mbs+1,PetscInt,&c->i);CHKERRQ(ierr); 2208 ierr = PetscMemcpy(c->i,a->i,(mbs+1)*sizeof(PetscInt));CHKERRQ(ierr); 2209 ierr = PetscLogObjectMemory((PetscObject)C,(mbs+1)*sizeof(PetscInt) + nz*(bs2*sizeof(MatScalar) + sizeof(PetscInt)));CHKERRQ(ierr); 2210 c->singlemalloc = PETSC_TRUE; 2211 c->free_a = PETSC_TRUE; 2212 c->free_ij = PETSC_TRUE; 2213 } 2214 if (mbs > 0) { 2215 if (cpvalues != MAT_SHARE_NONZERO_PATTERN) { 2216 ierr = PetscMemcpy(c->j,a->j,nz*sizeof(PetscInt));CHKERRQ(ierr); 2217 } 2218 if (cpvalues == MAT_COPY_VALUES) { 2219 ierr = PetscMemcpy(c->a,a->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2220 } else { 2221 ierr = PetscMemzero(c->a,bs2*nz*sizeof(MatScalar));CHKERRQ(ierr); 2222 } 2223 if (a->jshort) { 2224 /* cannot share jshort, it is reallocated in MatAssemblyEnd_SeqSBAIJ() */ 2225 /* if the parent matrix is reassembled, this child matrix will never notice */ 2226 ierr = PetscMalloc(nz*sizeof(unsigned short),&c->jshort);CHKERRQ(ierr); 2227 ierr = PetscLogObjectMemory((PetscObject)C,nz*sizeof(unsigned short));CHKERRQ(ierr); 2228 ierr = PetscMemcpy(c->jshort,a->jshort,nz*sizeof(unsigned short));CHKERRQ(ierr); 2229 2230 c->free_jshort = PETSC_TRUE; 2231 } 2232 } 2233 2234 c->roworiented = a->roworiented; 2235 c->nonew = a->nonew; 2236 2237 if (a->diag) { 2238 if (cpvalues == MAT_SHARE_NONZERO_PATTERN) { 2239 c->diag = a->diag; 2240 c->free_diag = PETSC_FALSE; 2241 } else { 2242 ierr = PetscMalloc(mbs*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 2243 ierr = PetscLogObjectMemory((PetscObject)C,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2244 for (i=0; i<mbs; i++) c->diag[i] = a->diag[i]; 2245 c->free_diag = PETSC_TRUE; 2246 } 2247 } 2248 c->nz = a->nz; 2249 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 2250 c->solve_work = 0; 2251 c->mult_work = 0; 2252 2253 *B = C; 2254 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 2255 PetscFunctionReturn(0); 2256 } 2257 2258 #undef __FUNCT__ 2259 #define __FUNCT__ "MatLoad_SeqSBAIJ" 2260 PetscErrorCode MatLoad_SeqSBAIJ(Mat newmat,PetscViewer viewer) 2261 { 2262 Mat_SeqSBAIJ *a; 2263 PetscErrorCode ierr; 2264 int fd; 2265 PetscMPIInt size; 2266 PetscInt i,nz,header[4],*rowlengths=0,M,N,bs=1; 2267 PetscInt *mask,mbs,*jj,j,rowcount,nzcount,k,*s_browlengths,maskcount; 2268 PetscInt kmax,jcount,block,idx,point,nzcountb,extra_rows,rows,cols; 2269 PetscInt *masked,nmask,tmp,bs2,ishift; 2270 PetscScalar *aa; 2271 MPI_Comm comm; 2272 2273 PetscFunctionBegin; 2274 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 2275 ierr = PetscOptionsGetInt(((PetscObject)newmat)->prefix,"-matload_block_size",&bs,NULL);CHKERRQ(ierr); 2276 bs2 = bs*bs; 2277 2278 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 2279 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"view must have one processor"); 2280 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 2281 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 2282 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not Mat object"); 2283 M = header[1]; N = header[2]; nz = header[3]; 2284 2285 if (header[3] < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format, cannot load as SeqSBAIJ"); 2286 2287 if (M != N) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Can only do square matrices"); 2288 2289 /* 2290 This code adds extra rows to make sure the number of rows is 2291 divisible by the blocksize 2292 */ 2293 mbs = M/bs; 2294 extra_rows = bs - M + bs*(mbs); 2295 if (extra_rows == bs) extra_rows = 0; 2296 else mbs++; 2297 if (extra_rows) { 2298 ierr = PetscInfo(viewer,"Padding loaded matrix to match blocksize\n");CHKERRQ(ierr); 2299 } 2300 2301 /* Set global sizes if not already set */ 2302 if (newmat->rmap->n < 0 && newmat->rmap->N < 0 && newmat->cmap->n < 0 && newmat->cmap->N < 0) { 2303 ierr = MatSetSizes(newmat,PETSC_DECIDE,PETSC_DECIDE,M+extra_rows,N+extra_rows);CHKERRQ(ierr); 2304 } else { /* Check if the matrix global sizes are correct */ 2305 ierr = MatGetSize(newmat,&rows,&cols);CHKERRQ(ierr); 2306 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); 2307 } 2308 2309 /* read in row lengths */ 2310 ierr = PetscMalloc((M+extra_rows)*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 2311 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 2312 for (i=0; i<extra_rows; i++) rowlengths[M+i] = 1; 2313 2314 /* read in column indices */ 2315 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscInt),&jj);CHKERRQ(ierr); 2316 ierr = PetscBinaryRead(fd,jj,nz,PETSC_INT);CHKERRQ(ierr); 2317 for (i=0; i<extra_rows; i++) jj[nz+i] = M+i; 2318 2319 /* loop over row lengths determining block row lengths */ 2320 ierr = PetscMalloc(mbs*sizeof(PetscInt),&s_browlengths);CHKERRQ(ierr); 2321 ierr = PetscMemzero(s_browlengths,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2322 ierr = PetscMalloc2(mbs,PetscInt,&mask,mbs,PetscInt,&masked);CHKERRQ(ierr); 2323 ierr = PetscMemzero(mask,mbs*sizeof(PetscInt));CHKERRQ(ierr); 2324 rowcount = 0; 2325 nzcount = 0; 2326 for (i=0; i<mbs; i++) { 2327 nmask = 0; 2328 for (j=0; j<bs; j++) { 2329 kmax = rowlengths[rowcount]; 2330 for (k=0; k<kmax; k++) { 2331 tmp = jj[nzcount++]/bs; /* block col. index */ 2332 if (!mask[tmp] && tmp >= i) {masked[nmask++] = tmp; mask[tmp] = 1;} 2333 } 2334 rowcount++; 2335 } 2336 s_browlengths[i] += nmask; 2337 2338 /* zero out the mask elements we set */ 2339 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2340 } 2341 2342 /* Do preallocation */ 2343 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(newmat,bs,0,s_browlengths);CHKERRQ(ierr); 2344 a = (Mat_SeqSBAIJ*)newmat->data; 2345 2346 /* set matrix "i" values */ 2347 a->i[0] = 0; 2348 for (i=1; i<= mbs; i++) { 2349 a->i[i] = a->i[i-1] + s_browlengths[i-1]; 2350 a->ilen[i-1] = s_browlengths[i-1]; 2351 } 2352 a->nz = a->i[mbs]; 2353 2354 /* read in nonzero values */ 2355 ierr = PetscMalloc((nz+extra_rows)*sizeof(PetscScalar),&aa);CHKERRQ(ierr); 2356 ierr = PetscBinaryRead(fd,aa,nz,PETSC_SCALAR);CHKERRQ(ierr); 2357 for (i=0; i<extra_rows; i++) aa[nz+i] = 1.0; 2358 2359 /* set "a" and "j" values into matrix */ 2360 nzcount = 0; jcount = 0; 2361 for (i=0; i<mbs; i++) { 2362 nzcountb = nzcount; 2363 nmask = 0; 2364 for (j=0; j<bs; j++) { 2365 kmax = rowlengths[i*bs+j]; 2366 for (k=0; k<kmax; k++) { 2367 tmp = jj[nzcount++]/bs; /* block col. index */ 2368 if (!mask[tmp] && tmp >= i) { masked[nmask++] = tmp; mask[tmp] = 1;} 2369 } 2370 } 2371 /* sort the masked values */ 2372 ierr = PetscSortInt(nmask,masked);CHKERRQ(ierr); 2373 2374 /* set "j" values into matrix */ 2375 maskcount = 1; 2376 for (j=0; j<nmask; j++) { 2377 a->j[jcount++] = masked[j]; 2378 mask[masked[j]] = maskcount++; 2379 } 2380 2381 /* set "a" values into matrix */ 2382 ishift = bs2*a->i[i]; 2383 for (j=0; j<bs; j++) { 2384 kmax = rowlengths[i*bs+j]; 2385 for (k=0; k<kmax; k++) { 2386 tmp = jj[nzcountb]/bs; /* block col. index */ 2387 if (tmp >= i) { 2388 block = mask[tmp] - 1; 2389 point = jj[nzcountb] - bs*tmp; 2390 idx = ishift + bs2*block + j + bs*point; 2391 a->a[idx] = aa[nzcountb]; 2392 } 2393 nzcountb++; 2394 } 2395 } 2396 /* zero out the mask elements we set */ 2397 for (j=0; j<nmask; j++) mask[masked[j]] = 0; 2398 } 2399 if (jcount != a->nz) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Bad binary matrix"); 2400 2401 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 2402 ierr = PetscFree(s_browlengths);CHKERRQ(ierr); 2403 ierr = PetscFree(aa);CHKERRQ(ierr); 2404 ierr = PetscFree(jj);CHKERRQ(ierr); 2405 ierr = PetscFree2(mask,masked);CHKERRQ(ierr); 2406 2407 ierr = MatAssemblyBegin(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2408 ierr = MatAssemblyEnd(newmat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2409 PetscFunctionReturn(0); 2410 } 2411 2412 #undef __FUNCT__ 2413 #define __FUNCT__ "MatCreateSeqSBAIJWithArrays" 2414 /*@ 2415 MatCreateSeqSBAIJWithArrays - Creates an sequential SBAIJ matrix using matrix elements 2416 (upper triangular entries in CSR format) provided by the user. 2417 2418 Collective on MPI_Comm 2419 2420 Input Parameters: 2421 + comm - must be an MPI communicator of size 1 2422 . bs - size of block 2423 . m - number of rows 2424 . n - number of columns 2425 . i - row indices 2426 . j - column indices 2427 - a - matrix values 2428 2429 Output Parameter: 2430 . mat - the matrix 2431 2432 Level: advanced 2433 2434 Notes: 2435 The i, j, and a arrays are not copied by this routine, the user must free these arrays 2436 once the matrix is destroyed 2437 2438 You cannot set new nonzero locations into this matrix, that will generate an error. 2439 2440 The i and j indices are 0 based 2441 2442 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 2443 it is the regular CSR format excluding the lower triangular elements. 2444 2445 .seealso: MatCreate(), MatCreateSBAIJ(), MatCreateSeqSBAIJ() 2446 2447 @*/ 2448 PetscErrorCode MatCreateSeqSBAIJWithArrays(MPI_Comm comm,PetscInt bs,PetscInt m,PetscInt n,PetscInt *i,PetscInt *j,PetscScalar *a,Mat *mat) 2449 { 2450 PetscErrorCode ierr; 2451 PetscInt ii; 2452 Mat_SeqSBAIJ *sbaij; 2453 2454 PetscFunctionBegin; 2455 if (bs != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"block size %D > 1 is not supported yet",bs); 2456 if (i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 2457 2458 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 2459 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 2460 ierr = MatSetType(*mat,MATSEQSBAIJ);CHKERRQ(ierr); 2461 ierr = MatSeqSBAIJSetPreallocation_SeqSBAIJ(*mat,bs,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 2462 sbaij = (Mat_SeqSBAIJ*)(*mat)->data; 2463 ierr = PetscMalloc2(m,PetscInt,&sbaij->imax,m,PetscInt,&sbaij->ilen);CHKERRQ(ierr); 2464 ierr = PetscLogObjectMemory((PetscObject)*mat,2*m*sizeof(PetscInt));CHKERRQ(ierr); 2465 2466 sbaij->i = i; 2467 sbaij->j = j; 2468 sbaij->a = a; 2469 2470 sbaij->singlemalloc = PETSC_FALSE; 2471 sbaij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 2472 sbaij->free_a = PETSC_FALSE; 2473 sbaij->free_ij = PETSC_FALSE; 2474 2475 for (ii=0; ii<m; ii++) { 2476 sbaij->ilen[ii] = sbaij->imax[ii] = i[ii+1] - i[ii]; 2477 #if defined(PETSC_USE_DEBUG) 2478 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]); 2479 #endif 2480 } 2481 #if defined(PETSC_USE_DEBUG) 2482 for (ii=0; ii<sbaij->i[m]; ii++) { 2483 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 2484 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]); 2485 } 2486 #endif 2487 2488 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2489 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2490 PetscFunctionReturn(0); 2491 } 2492 2493 2494 2495 2496 2497