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