1 /* 2 Defines the basic matrix operations for the AIJ (compressed row) 3 matrix storage format. 4 */ 5 6 #include <../src/mat/impls/aij/seq/aij.h> /*I "petscmat.h" I*/ 7 #include <petscblaslapack.h> 8 #include <petscbt.h> 9 #include <petsc/private/kernels/blocktranspose.h> 10 11 PetscErrorCode MatSeqAIJSetTypeFromOptions(Mat A) 12 { 13 PetscErrorCode ierr; 14 PetscBool flg; 15 char type[256]; 16 17 PetscFunctionBegin; 18 ierr = PetscObjectOptionsBegin((PetscObject)A);PetscCall(ierr); 19 PetscCall(PetscOptionsFList("-mat_seqaij_type","Matrix SeqAIJ type","MatSeqAIJSetType",MatSeqAIJList,"seqaij",type,256,&flg)); 20 if (flg) PetscCall(MatSeqAIJSetType(A,type)); 21 ierr = PetscOptionsEnd();PetscCall(ierr); 22 PetscFunctionReturn(0); 23 } 24 25 PetscErrorCode MatGetColumnReductions_SeqAIJ(Mat A,PetscInt type,PetscReal *reductions) 26 { 27 PetscInt i,m,n; 28 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 29 30 PetscFunctionBegin; 31 PetscCall(MatGetSize(A,&m,&n)); 32 PetscCall(PetscArrayzero(reductions,n)); 33 if (type == NORM_2) { 34 for (i=0; i<aij->i[m]; i++) { 35 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]*aij->a[i]); 36 } 37 } else if (type == NORM_1) { 38 for (i=0; i<aij->i[m]; i++) { 39 reductions[aij->j[i]] += PetscAbsScalar(aij->a[i]); 40 } 41 } else if (type == NORM_INFINITY) { 42 for (i=0; i<aij->i[m]; i++) { 43 reductions[aij->j[i]] = PetscMax(PetscAbsScalar(aij->a[i]),reductions[aij->j[i]]); 44 } 45 } else if (type == REDUCTION_SUM_REALPART || type == REDUCTION_MEAN_REALPART) { 46 for (i=0; i<aij->i[m]; i++) { 47 reductions[aij->j[i]] += PetscRealPart(aij->a[i]); 48 } 49 } else if (type == REDUCTION_SUM_IMAGINARYPART || type == REDUCTION_MEAN_IMAGINARYPART) { 50 for (i=0; i<aij->i[m]; i++) { 51 reductions[aij->j[i]] += PetscImaginaryPart(aij->a[i]); 52 } 53 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Unknown reduction type"); 54 55 if (type == NORM_2) { 56 for (i=0; i<n; i++) reductions[i] = PetscSqrtReal(reductions[i]); 57 } else if (type == REDUCTION_MEAN_REALPART || type == REDUCTION_MEAN_IMAGINARYPART) { 58 for (i=0; i<n; i++) reductions[i] /= m; 59 } 60 PetscFunctionReturn(0); 61 } 62 63 PetscErrorCode MatFindOffBlockDiagonalEntries_SeqAIJ(Mat A,IS *is) 64 { 65 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 66 PetscInt i,m=A->rmap->n,cnt = 0, bs = A->rmap->bs; 67 const PetscInt *jj = a->j,*ii = a->i; 68 PetscInt *rows; 69 70 PetscFunctionBegin; 71 for (i=0; i<m; i++) { 72 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 73 cnt++; 74 } 75 } 76 PetscCall(PetscMalloc1(cnt,&rows)); 77 cnt = 0; 78 for (i=0; i<m; i++) { 79 if ((ii[i] != ii[i+1]) && ((jj[ii[i]] < bs*(i/bs)) || (jj[ii[i+1]-1] > bs*((i+bs)/bs)-1))) { 80 rows[cnt] = i; 81 cnt++; 82 } 83 } 84 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,is)); 85 PetscFunctionReturn(0); 86 } 87 88 PetscErrorCode MatFindZeroDiagonals_SeqAIJ_Private(Mat A,PetscInt *nrows,PetscInt **zrows) 89 { 90 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 91 const MatScalar *aa; 92 PetscInt i,m=A->rmap->n,cnt = 0; 93 const PetscInt *ii = a->i,*jj = a->j,*diag; 94 PetscInt *rows; 95 96 PetscFunctionBegin; 97 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 98 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 99 diag = a->diag; 100 for (i=0; i<m; i++) { 101 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 102 cnt++; 103 } 104 } 105 PetscCall(PetscMalloc1(cnt,&rows)); 106 cnt = 0; 107 for (i=0; i<m; i++) { 108 if ((diag[i] >= ii[i+1]) || (jj[diag[i]] != i) || (aa[diag[i]] == 0.0)) { 109 rows[cnt++] = i; 110 } 111 } 112 *nrows = cnt; 113 *zrows = rows; 114 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 115 PetscFunctionReturn(0); 116 } 117 118 PetscErrorCode MatFindZeroDiagonals_SeqAIJ(Mat A,IS *zrows) 119 { 120 PetscInt nrows,*rows; 121 122 PetscFunctionBegin; 123 *zrows = NULL; 124 PetscCall(MatFindZeroDiagonals_SeqAIJ_Private(A,&nrows,&rows)); 125 PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)A),nrows,rows,PETSC_OWN_POINTER,zrows)); 126 PetscFunctionReturn(0); 127 } 128 129 PetscErrorCode MatFindNonzeroRows_SeqAIJ(Mat A,IS *keptrows) 130 { 131 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 132 const MatScalar *aa; 133 PetscInt m=A->rmap->n,cnt = 0; 134 const PetscInt *ii; 135 PetscInt n,i,j,*rows; 136 137 PetscFunctionBegin; 138 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 139 *keptrows = NULL; 140 ii = a->i; 141 for (i=0; i<m; i++) { 142 n = ii[i+1] - ii[i]; 143 if (!n) { 144 cnt++; 145 goto ok1; 146 } 147 for (j=ii[i]; j<ii[i+1]; j++) { 148 if (aa[j] != 0.0) goto ok1; 149 } 150 cnt++; 151 ok1:; 152 } 153 if (!cnt) { 154 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 155 PetscFunctionReturn(0); 156 } 157 PetscCall(PetscMalloc1(A->rmap->n-cnt,&rows)); 158 cnt = 0; 159 for (i=0; i<m; i++) { 160 n = ii[i+1] - ii[i]; 161 if (!n) continue; 162 for (j=ii[i]; j<ii[i+1]; j++) { 163 if (aa[j] != 0.0) { 164 rows[cnt++] = i; 165 break; 166 } 167 } 168 } 169 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 170 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,cnt,rows,PETSC_OWN_POINTER,keptrows)); 171 PetscFunctionReturn(0); 172 } 173 174 PetscErrorCode MatDiagonalSet_SeqAIJ(Mat Y,Vec D,InsertMode is) 175 { 176 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) Y->data; 177 PetscInt i,m = Y->rmap->n; 178 const PetscInt *diag; 179 MatScalar *aa; 180 const PetscScalar *v; 181 PetscBool missing; 182 183 PetscFunctionBegin; 184 if (Y->assembled) { 185 PetscCall(MatMissingDiagonal_SeqAIJ(Y,&missing,NULL)); 186 if (!missing) { 187 diag = aij->diag; 188 PetscCall(VecGetArrayRead(D,&v)); 189 PetscCall(MatSeqAIJGetArray(Y,&aa)); 190 if (is == INSERT_VALUES) { 191 for (i=0; i<m; i++) { 192 aa[diag[i]] = v[i]; 193 } 194 } else { 195 for (i=0; i<m; i++) { 196 aa[diag[i]] += v[i]; 197 } 198 } 199 PetscCall(MatSeqAIJRestoreArray(Y,&aa)); 200 PetscCall(VecRestoreArrayRead(D,&v)); 201 PetscFunctionReturn(0); 202 } 203 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 204 } 205 PetscCall(MatDiagonalSet_Default(Y,D,is)); 206 PetscFunctionReturn(0); 207 } 208 209 PetscErrorCode MatGetRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *m,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 210 { 211 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 212 PetscInt i,ishift; 213 214 PetscFunctionBegin; 215 *m = A->rmap->n; 216 if (!ia) PetscFunctionReturn(0); 217 ishift = 0; 218 if (symmetric && !A->structurally_symmetric) { 219 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,ishift,oshift,(PetscInt**)ia,(PetscInt**)ja)); 220 } else if (oshift == 1) { 221 PetscInt *tia; 222 PetscInt nz = a->i[A->rmap->n]; 223 /* malloc space and add 1 to i and j indices */ 224 PetscCall(PetscMalloc1(A->rmap->n+1,&tia)); 225 for (i=0; i<A->rmap->n+1; i++) tia[i] = a->i[i] + 1; 226 *ia = tia; 227 if (ja) { 228 PetscInt *tja; 229 PetscCall(PetscMalloc1(nz+1,&tja)); 230 for (i=0; i<nz; i++) tja[i] = a->j[i] + 1; 231 *ja = tja; 232 } 233 } else { 234 *ia = a->i; 235 if (ja) *ja = a->j; 236 } 237 PetscFunctionReturn(0); 238 } 239 240 PetscErrorCode MatRestoreRowIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 241 { 242 PetscFunctionBegin; 243 if (!ia) PetscFunctionReturn(0); 244 if ((symmetric && !A->structurally_symmetric) || oshift == 1) { 245 PetscCall(PetscFree(*ia)); 246 if (ja) PetscCall(PetscFree(*ja)); 247 } 248 PetscFunctionReturn(0); 249 } 250 251 PetscErrorCode MatGetColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 252 { 253 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 254 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 255 PetscInt nz = a->i[m],row,*jj,mr,col; 256 257 PetscFunctionBegin; 258 *nn = n; 259 if (!ia) PetscFunctionReturn(0); 260 if (symmetric) { 261 PetscCall(MatToSymmetricIJ_SeqAIJ(A->rmap->n,a->i,a->j,PETSC_TRUE,0,oshift,(PetscInt**)ia,(PetscInt**)ja)); 262 } else { 263 PetscCall(PetscCalloc1(n,&collengths)); 264 PetscCall(PetscMalloc1(n+1,&cia)); 265 PetscCall(PetscMalloc1(nz,&cja)); 266 jj = a->j; 267 for (i=0; i<nz; i++) { 268 collengths[jj[i]]++; 269 } 270 cia[0] = oshift; 271 for (i=0; i<n; i++) { 272 cia[i+1] = cia[i] + collengths[i]; 273 } 274 PetscCall(PetscArrayzero(collengths,n)); 275 jj = a->j; 276 for (row=0; row<m; row++) { 277 mr = a->i[row+1] - a->i[row]; 278 for (i=0; i<mr; i++) { 279 col = *jj++; 280 281 cja[cia[col] + collengths[col]++ - oshift] = row + oshift; 282 } 283 } 284 PetscCall(PetscFree(collengths)); 285 *ia = cia; *ja = cja; 286 } 287 PetscFunctionReturn(0); 288 } 289 290 PetscErrorCode MatRestoreColumnIJ_SeqAIJ(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscBool *done) 291 { 292 PetscFunctionBegin; 293 if (!ia) PetscFunctionReturn(0); 294 295 PetscCall(PetscFree(*ia)); 296 PetscCall(PetscFree(*ja)); 297 PetscFunctionReturn(0); 298 } 299 300 /* 301 MatGetColumnIJ_SeqAIJ_Color() and MatRestoreColumnIJ_SeqAIJ_Color() are customized from 302 MatGetColumnIJ_SeqAIJ() and MatRestoreColumnIJ_SeqAIJ() by adding an output 303 spidx[], index of a->a, to be used in MatTransposeColoringCreate_SeqAIJ() and MatFDColoringCreate_SeqXAIJ() 304 */ 305 PetscErrorCode MatGetColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *nn,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 306 { 307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 308 PetscInt i,*collengths,*cia,*cja,n = A->cmap->n,m = A->rmap->n; 309 PetscInt nz = a->i[m],row,mr,col,tmp; 310 PetscInt *cspidx; 311 const PetscInt *jj; 312 313 PetscFunctionBegin; 314 *nn = n; 315 if (!ia) PetscFunctionReturn(0); 316 317 PetscCall(PetscCalloc1(n,&collengths)); 318 PetscCall(PetscMalloc1(n+1,&cia)); 319 PetscCall(PetscMalloc1(nz,&cja)); 320 PetscCall(PetscMalloc1(nz,&cspidx)); 321 jj = a->j; 322 for (i=0; i<nz; i++) { 323 collengths[jj[i]]++; 324 } 325 cia[0] = oshift; 326 for (i=0; i<n; i++) { 327 cia[i+1] = cia[i] + collengths[i]; 328 } 329 PetscCall(PetscArrayzero(collengths,n)); 330 jj = a->j; 331 for (row=0; row<m; row++) { 332 mr = a->i[row+1] - a->i[row]; 333 for (i=0; i<mr; i++) { 334 col = *jj++; 335 tmp = cia[col] + collengths[col]++ - oshift; 336 cspidx[tmp] = a->i[row] + i; /* index of a->j */ 337 cja[tmp] = row + oshift; 338 } 339 } 340 PetscCall(PetscFree(collengths)); 341 *ia = cia; 342 *ja = cja; 343 *spidx = cspidx; 344 PetscFunctionReturn(0); 345 } 346 347 PetscErrorCode MatRestoreColumnIJ_SeqAIJ_Color(Mat A,PetscInt oshift,PetscBool symmetric,PetscBool inodecompressed,PetscInt *n,const PetscInt *ia[],const PetscInt *ja[],PetscInt *spidx[],PetscBool *done) 348 { 349 PetscFunctionBegin; 350 PetscCall(MatRestoreColumnIJ_SeqAIJ(A,oshift,symmetric,inodecompressed,n,ia,ja,done)); 351 PetscCall(PetscFree(*spidx)); 352 PetscFunctionReturn(0); 353 } 354 355 PetscErrorCode MatSetValuesRow_SeqAIJ(Mat A,PetscInt row,const PetscScalar v[]) 356 { 357 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 358 PetscInt *ai = a->i; 359 PetscScalar *aa; 360 361 PetscFunctionBegin; 362 PetscCall(MatSeqAIJGetArray(A,&aa)); 363 PetscCall(PetscArraycpy(aa+ai[row],v,ai[row+1]-ai[row])); 364 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 365 PetscFunctionReturn(0); 366 } 367 368 /* 369 MatSeqAIJSetValuesLocalFast - An optimized version of MatSetValuesLocal() for SeqAIJ matrices with several assumptions 370 371 - a single row of values is set with each call 372 - no row or column indices are negative or (in error) larger than the number of rows or columns 373 - the values are always added to the matrix, not set 374 - no new locations are introduced in the nonzero structure of the matrix 375 376 This does NOT assume the global column indices are sorted 377 378 */ 379 380 #include <petsc/private/isimpl.h> 381 PetscErrorCode MatSeqAIJSetValuesLocalFast(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 382 { 383 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 384 PetscInt low,high,t,row,nrow,i,col,l; 385 const PetscInt *rp,*ai = a->i,*ailen = a->ilen,*aj = a->j; 386 PetscInt lastcol = -1; 387 MatScalar *ap,value,*aa; 388 const PetscInt *ridx = A->rmap->mapping->indices,*cidx = A->cmap->mapping->indices; 389 390 PetscFunctionBegin; 391 PetscCall(MatSeqAIJGetArray(A,&aa)); 392 row = ridx[im[0]]; 393 rp = aj + ai[row]; 394 ap = aa + ai[row]; 395 nrow = ailen[row]; 396 low = 0; 397 high = nrow; 398 for (l=0; l<n; l++) { /* loop over added columns */ 399 col = cidx[in[l]]; 400 value = v[l]; 401 402 if (col <= lastcol) low = 0; 403 else high = nrow; 404 lastcol = col; 405 while (high-low > 5) { 406 t = (low+high)/2; 407 if (rp[t] > col) high = t; 408 else low = t; 409 } 410 for (i=low; i<high; i++) { 411 if (rp[i] == col) { 412 ap[i] += value; 413 low = i + 1; 414 break; 415 } 416 } 417 } 418 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 419 return 0; 420 } 421 422 PetscErrorCode MatSetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 423 { 424 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 425 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 426 PetscInt *imax = a->imax,*ai = a->i,*ailen = a->ilen; 427 PetscInt *aj = a->j,nonew = a->nonew,lastcol = -1; 428 MatScalar *ap=NULL,value=0.0,*aa; 429 PetscBool ignorezeroentries = a->ignorezeroentries; 430 PetscBool roworiented = a->roworiented; 431 432 PetscFunctionBegin; 433 PetscCall(MatSeqAIJGetArray(A,&aa)); 434 for (k=0; k<m; k++) { /* loop over added rows */ 435 row = im[k]; 436 if (row < 0) continue; 437 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 438 rp = aj + ai[row]; 439 if (!A->structure_only) ap = aa + ai[row]; 440 rmax = imax[row]; nrow = ailen[row]; 441 low = 0; 442 high = nrow; 443 for (l=0; l<n; l++) { /* loop over added columns */ 444 if (in[l] < 0) continue; 445 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 446 col = in[l]; 447 if (v && !A->structure_only) value = roworiented ? v[l + k*n] : v[k + l*m]; 448 if (!A->structure_only && value == 0.0 && ignorezeroentries && is == ADD_VALUES && row != col) continue; 449 450 if (col <= lastcol) low = 0; 451 else high = nrow; 452 lastcol = col; 453 while (high-low > 5) { 454 t = (low+high)/2; 455 if (rp[t] > col) high = t; 456 else low = t; 457 } 458 for (i=low; i<high; i++) { 459 if (rp[i] > col) break; 460 if (rp[i] == col) { 461 if (!A->structure_only) { 462 if (is == ADD_VALUES) { 463 ap[i] += value; 464 (void)PetscLogFlops(1.0); 465 } 466 else ap[i] = value; 467 } 468 low = i + 1; 469 goto noinsert; 470 } 471 } 472 if (value == 0.0 && ignorezeroentries && row != col) goto noinsert; 473 if (nonew == 1) goto noinsert; 474 PetscCheckFalse(nonew == -1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero at (%" PetscInt_FMT ",%" PetscInt_FMT ") in the matrix",row,col); 475 if (A->structure_only) { 476 MatSeqXAIJReallocateAIJ_structure_only(A,A->rmap->n,1,nrow,row,col,rmax,ai,aj,rp,imax,nonew,MatScalar); 477 } else { 478 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 479 } 480 N = nrow++ - 1; a->nz++; high++; 481 /* shift up all the later entries in this row */ 482 PetscCall(PetscArraymove(rp+i+1,rp+i,N-i+1)); 483 rp[i] = col; 484 if (!A->structure_only) { 485 PetscCall(PetscArraymove(ap+i+1,ap+i,N-i+1)); 486 ap[i] = value; 487 } 488 low = i + 1; 489 A->nonzerostate++; 490 noinsert:; 491 } 492 ailen[row] = nrow; 493 } 494 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 495 PetscFunctionReturn(0); 496 } 497 498 PetscErrorCode MatSetValues_SeqAIJ_SortedFullNoPreallocation(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 499 { 500 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 501 PetscInt *rp,k,row; 502 PetscInt *ai = a->i; 503 PetscInt *aj = a->j; 504 MatScalar *aa,*ap; 505 506 PetscFunctionBegin; 507 PetscCheck(!A->was_assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot call on assembled matrix."); 508 PetscCheckFalse(m*n+a->nz > a->maxnz,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Number of entries in matrix will be larger than maximum nonzeros allocated for %" PetscInt_FMT " in MatSeqAIJSetTotalPreallocation()",a->maxnz); 509 510 PetscCall(MatSeqAIJGetArray(A,&aa)); 511 for (k=0; k<m; k++) { /* loop over added rows */ 512 row = im[k]; 513 rp = aj + ai[row]; 514 ap = aa + ai[row]; 515 516 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 517 if (!A->structure_only) { 518 if (v) { 519 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 520 v += n; 521 } else { 522 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 523 } 524 } 525 a->ilen[row] = n; 526 a->imax[row] = n; 527 a->i[row+1] = a->i[row]+n; 528 a->nz += n; 529 } 530 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 531 PetscFunctionReturn(0); 532 } 533 534 /*@ 535 MatSeqAIJSetTotalPreallocation - Sets an upper bound on the total number of expected nonzeros in the matrix. 536 537 Input Parameters: 538 + A - the SeqAIJ matrix 539 - nztotal - bound on the number of nonzeros 540 541 Level: advanced 542 543 Notes: 544 This can be called if you will be provided the matrix row by row (from row zero) with sorted column indices for each row. 545 Simply call MatSetValues() after this call to provide the matrix entries in the usual manner. This matrix may be used 546 as always with multiple matrix assemblies. 547 548 .seealso: MatSetOption(), MAT_SORTED_FULL, MatSetValues(), MatSeqAIJSetPreallocation() 549 @*/ 550 551 PetscErrorCode MatSeqAIJSetTotalPreallocation(Mat A,PetscInt nztotal) 552 { 553 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 554 555 PetscFunctionBegin; 556 PetscCall(PetscLayoutSetUp(A->rmap)); 557 PetscCall(PetscLayoutSetUp(A->cmap)); 558 a->maxnz = nztotal; 559 if (!a->imax) { 560 PetscCall(PetscMalloc1(A->rmap->n,&a->imax)); 561 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 562 } 563 if (!a->ilen) { 564 PetscCall(PetscMalloc1(A->rmap->n,&a->ilen)); 565 PetscCall(PetscLogObjectMemory((PetscObject)A,A->rmap->n*sizeof(PetscInt))); 566 } else { 567 PetscCall(PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt))); 568 } 569 570 /* allocate the matrix space */ 571 if (A->structure_only) { 572 PetscCall(PetscMalloc1(nztotal,&a->j)); 573 PetscCall(PetscMalloc1(A->rmap->n+1,&a->i)); 574 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*sizeof(PetscInt))); 575 } else { 576 PetscCall(PetscMalloc3(nztotal,&a->a,nztotal,&a->j,A->rmap->n+1,&a->i)); 577 PetscCall(PetscLogObjectMemory((PetscObject)A,(A->rmap->n+1)*sizeof(PetscInt)+nztotal*(sizeof(PetscScalar)+sizeof(PetscInt)))); 578 } 579 a->i[0] = 0; 580 if (A->structure_only) { 581 a->singlemalloc = PETSC_FALSE; 582 a->free_a = PETSC_FALSE; 583 } else { 584 a->singlemalloc = PETSC_TRUE; 585 a->free_a = PETSC_TRUE; 586 } 587 a->free_ij = PETSC_TRUE; 588 A->ops->setvalues = MatSetValues_SeqAIJ_SortedFullNoPreallocation; 589 A->preallocated = PETSC_TRUE; 590 PetscFunctionReturn(0); 591 } 592 593 PetscErrorCode MatSetValues_SeqAIJ_SortedFull(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],const PetscScalar v[],InsertMode is) 594 { 595 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 596 PetscInt *rp,k,row; 597 PetscInt *ai = a->i,*ailen = a->ilen; 598 PetscInt *aj = a->j; 599 MatScalar *aa,*ap; 600 601 PetscFunctionBegin; 602 PetscCall(MatSeqAIJGetArray(A,&aa)); 603 for (k=0; k<m; k++) { /* loop over added rows */ 604 row = im[k]; 605 PetscCheck(n <= a->imax[row],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Preallocation for row %" PetscInt_FMT " does not match number of columns provided",n); 606 rp = aj + ai[row]; 607 ap = aa + ai[row]; 608 if (!A->was_assembled) { 609 PetscCall(PetscMemcpy(rp,in,n*sizeof(PetscInt))); 610 } 611 if (!A->structure_only) { 612 if (v) { 613 PetscCall(PetscMemcpy(ap,v,n*sizeof(PetscScalar))); 614 v += n; 615 } else { 616 PetscCall(PetscMemzero(ap,n*sizeof(PetscScalar))); 617 } 618 } 619 ailen[row] = n; 620 a->nz += n; 621 } 622 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 623 PetscFunctionReturn(0); 624 } 625 626 PetscErrorCode MatGetValues_SeqAIJ(Mat A,PetscInt m,const PetscInt im[],PetscInt n,const PetscInt in[],PetscScalar v[]) 627 { 628 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 629 PetscInt *rp,k,low,high,t,row,nrow,i,col,l,*aj = a->j; 630 PetscInt *ai = a->i,*ailen = a->ilen; 631 MatScalar *ap,*aa; 632 633 PetscFunctionBegin; 634 PetscCall(MatSeqAIJGetArray(A,&aa)); 635 for (k=0; k<m; k++) { /* loop over rows */ 636 row = im[k]; 637 if (row < 0) {v += n; continue;} /* negative row */ 638 PetscCheck(row < A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row too large: row %" PetscInt_FMT " max %" PetscInt_FMT,row,A->rmap->n-1); 639 rp = aj + ai[row]; ap = aa + ai[row]; 640 nrow = ailen[row]; 641 for (l=0; l<n; l++) { /* loop over columns */ 642 if (in[l] < 0) {v++; continue;} /* negative column */ 643 PetscCheck(in[l] < A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column too large: col %" PetscInt_FMT " max %" PetscInt_FMT,in[l],A->cmap->n-1); 644 col = in[l]; 645 high = nrow; low = 0; /* assume unsorted */ 646 while (high-low > 5) { 647 t = (low+high)/2; 648 if (rp[t] > col) high = t; 649 else low = t; 650 } 651 for (i=low; i<high; i++) { 652 if (rp[i] > col) break; 653 if (rp[i] == col) { 654 *v++ = ap[i]; 655 goto finished; 656 } 657 } 658 *v++ = 0.0; 659 finished:; 660 } 661 } 662 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 663 PetscFunctionReturn(0); 664 } 665 666 PetscErrorCode MatView_SeqAIJ_Binary(Mat mat,PetscViewer viewer) 667 { 668 Mat_SeqAIJ *A = (Mat_SeqAIJ*)mat->data; 669 const PetscScalar *av; 670 PetscInt header[4],M,N,m,nz,i; 671 PetscInt *rowlens; 672 673 PetscFunctionBegin; 674 PetscCall(PetscViewerSetUp(viewer)); 675 676 M = mat->rmap->N; 677 N = mat->cmap->N; 678 m = mat->rmap->n; 679 nz = A->nz; 680 681 /* write matrix header */ 682 header[0] = MAT_FILE_CLASSID; 683 header[1] = M; header[2] = N; header[3] = nz; 684 PetscCall(PetscViewerBinaryWrite(viewer,header,4,PETSC_INT)); 685 686 /* fill in and store row lengths */ 687 PetscCall(PetscMalloc1(m,&rowlens)); 688 for (i=0; i<m; i++) rowlens[i] = A->i[i+1] - A->i[i]; 689 PetscCall(PetscViewerBinaryWrite(viewer,rowlens,m,PETSC_INT)); 690 PetscCall(PetscFree(rowlens)); 691 /* store column indices */ 692 PetscCall(PetscViewerBinaryWrite(viewer,A->j,nz,PETSC_INT)); 693 /* store nonzero values */ 694 PetscCall(MatSeqAIJGetArrayRead(mat,&av)); 695 PetscCall(PetscViewerBinaryWrite(viewer,av,nz,PETSC_SCALAR)); 696 PetscCall(MatSeqAIJRestoreArrayRead(mat,&av)); 697 698 /* write block size option to the viewer's .info file */ 699 PetscCall(MatView_Binary_BlockSizes(mat,viewer)); 700 PetscFunctionReturn(0); 701 } 702 703 static PetscErrorCode MatView_SeqAIJ_ASCII_structonly(Mat A,PetscViewer viewer) 704 { 705 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 706 PetscInt i,k,m=A->rmap->N; 707 708 PetscFunctionBegin; 709 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 710 for (i=0; i<m; i++) { 711 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 712 for (k=a->i[i]; k<a->i[i+1]; k++) { 713 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ") ",a->j[k])); 714 } 715 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 716 } 717 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 718 PetscFunctionReturn(0); 719 } 720 721 extern PetscErrorCode MatSeqAIJFactorInfo_Matlab(Mat,PetscViewer); 722 723 PetscErrorCode MatView_SeqAIJ_ASCII(Mat A,PetscViewer viewer) 724 { 725 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 726 const PetscScalar *av; 727 PetscInt i,j,m = A->rmap->n; 728 const char *name; 729 PetscViewerFormat format; 730 731 PetscFunctionBegin; 732 if (A->structure_only) { 733 PetscCall(MatView_SeqAIJ_ASCII_structonly(A,viewer)); 734 PetscFunctionReturn(0); 735 } 736 737 PetscCall(PetscViewerGetFormat(viewer,&format)); 738 if (format == PETSC_VIEWER_ASCII_FACTOR_INFO || format == PETSC_VIEWER_ASCII_INFO || format == PETSC_VIEWER_ASCII_INFO_DETAIL) PetscFunctionReturn(0); 739 740 /* trigger copy to CPU if needed */ 741 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 742 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 743 if (format == PETSC_VIEWER_ASCII_MATLAB) { 744 PetscInt nofinalvalue = 0; 745 if (m && ((a->i[m] == a->i[m-1]) || (a->j[a->nz-1] != A->cmap->n-1))) { 746 /* Need a dummy value to ensure the dimension of the matrix. */ 747 nofinalvalue = 1; 748 } 749 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 750 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Size = %" PetscInt_FMT " %" PetscInt_FMT " \n",m,A->cmap->n)); 751 PetscCall(PetscViewerASCIIPrintf(viewer,"%% Nonzeros = %" PetscInt_FMT " \n",a->nz)); 752 #if defined(PETSC_USE_COMPLEX) 753 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",4);\n",a->nz+nofinalvalue)); 754 #else 755 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = zeros(%" PetscInt_FMT ",3);\n",a->nz+nofinalvalue)); 756 #endif 757 PetscCall(PetscViewerASCIIPrintf(viewer,"zzz = [\n")); 758 759 for (i=0; i<m; i++) { 760 for (j=a->i[i]; j<a->i[i+1]; j++) { 761 #if defined(PETSC_USE_COMPLEX) 762 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",i+1,a->j[j]+1,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 763 #else 764 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",i+1,a->j[j]+1,(double)a->a[j])); 765 #endif 766 } 767 } 768 if (nofinalvalue) { 769 #if defined(PETSC_USE_COMPLEX) 770 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e %18.16e\n",m,A->cmap->n,0.,0.)); 771 #else 772 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %18.16e\n",m,A->cmap->n,0.0)); 773 #endif 774 } 775 PetscCall(PetscObjectGetName((PetscObject)A,&name)); 776 PetscCall(PetscViewerASCIIPrintf(viewer,"];\n %s = spconvert(zzz);\n",name)); 777 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 778 } else if (format == PETSC_VIEWER_ASCII_COMMON) { 779 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 780 for (i=0; i<m; i++) { 781 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 782 for (j=a->i[i]; j<a->i[i+1]; j++) { 783 #if defined(PETSC_USE_COMPLEX) 784 if (PetscImaginaryPart(a->a[j]) > 0.0 && PetscRealPart(a->a[j]) != 0.0) { 785 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 786 } else if (PetscImaginaryPart(a->a[j]) < 0.0 && PetscRealPart(a->a[j]) != 0.0) { 787 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 788 } else if (PetscRealPart(a->a[j]) != 0.0) { 789 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 790 } 791 #else 792 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 793 #endif 794 } 795 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 796 } 797 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 798 } else if (format == PETSC_VIEWER_ASCII_SYMMODU) { 799 PetscInt nzd=0,fshift=1,*sptr; 800 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 801 PetscCall(PetscMalloc1(m+1,&sptr)); 802 for (i=0; i<m; i++) { 803 sptr[i] = nzd+1; 804 for (j=a->i[i]; j<a->i[i+1]; j++) { 805 if (a->j[j] >= i) { 806 #if defined(PETSC_USE_COMPLEX) 807 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) nzd++; 808 #else 809 if (a->a[j] != 0.0) nzd++; 810 #endif 811 } 812 } 813 } 814 sptr[m] = nzd+1; 815 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n\n",m,nzd)); 816 for (i=0; i<m+1; i+=6) { 817 if (i+4<m) { 818 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4],sptr[i+5])); 819 } else if (i+3<m) { 820 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3],sptr[i+4])); 821 } else if (i+2<m) { 822 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2],sptr[i+3])); 823 } else if (i+1<m) { 824 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1],sptr[i+2])); 825 } else if (i<m) { 826 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " %" PetscInt_FMT "\n",sptr[i],sptr[i+1])); 827 } else { 828 PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT "\n",sptr[i])); 829 } 830 } 831 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 832 PetscCall(PetscFree(sptr)); 833 for (i=0; i<m; i++) { 834 for (j=a->i[i]; j<a->i[i+1]; j++) { 835 if (a->j[j] >= i) PetscCall(PetscViewerASCIIPrintf(viewer," %" PetscInt_FMT " ",a->j[j]+fshift)); 836 } 837 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 838 } 839 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 840 for (i=0; i<m; i++) { 841 for (j=a->i[i]; j<a->i[i+1]; j++) { 842 if (a->j[j] >= i) { 843 #if defined(PETSC_USE_COMPLEX) 844 if (PetscImaginaryPart(a->a[j]) != 0.0 || PetscRealPart(a->a[j]) != 0.0) { 845 PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e %18.16e ",(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 846 } 847 #else 848 if (a->a[j] != 0.0) PetscCall(PetscViewerASCIIPrintf(viewer," %18.16e ",(double)a->a[j])); 849 #endif 850 } 851 } 852 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 853 } 854 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 855 } else if (format == PETSC_VIEWER_ASCII_DENSE) { 856 PetscInt cnt = 0,jcnt; 857 PetscScalar value; 858 #if defined(PETSC_USE_COMPLEX) 859 PetscBool realonly = PETSC_TRUE; 860 861 for (i=0; i<a->i[m]; i++) { 862 if (PetscImaginaryPart(a->a[i]) != 0.0) { 863 realonly = PETSC_FALSE; 864 break; 865 } 866 } 867 #endif 868 869 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 870 for (i=0; i<m; i++) { 871 jcnt = 0; 872 for (j=0; j<A->cmap->n; j++) { 873 if (jcnt < a->i[i+1]-a->i[i] && j == a->j[cnt]) { 874 value = a->a[cnt++]; 875 jcnt++; 876 } else { 877 value = 0.0; 878 } 879 #if defined(PETSC_USE_COMPLEX) 880 if (realonly) { 881 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)PetscRealPart(value))); 882 } else { 883 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e+%7.5e i ",(double)PetscRealPart(value),(double)PetscImaginaryPart(value))); 884 } 885 #else 886 PetscCall(PetscViewerASCIIPrintf(viewer," %7.5e ",(double)value)); 887 #endif 888 } 889 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 890 } 891 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 892 } else if (format == PETSC_VIEWER_ASCII_MATRIXMARKET) { 893 PetscInt fshift=1; 894 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 895 #if defined(PETSC_USE_COMPLEX) 896 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate complex general\n")); 897 #else 898 PetscCall(PetscViewerASCIIPrintf(viewer,"%%%%MatrixMarket matrix coordinate real general\n")); 899 #endif 900 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %" PetscInt_FMT "\n", m, A->cmap->n, a->nz)); 901 for (i=0; i<m; i++) { 902 for (j=a->i[i]; j<a->i[i+1]; j++) { 903 #if defined(PETSC_USE_COMPLEX) 904 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g %g\n", i+fshift,a->j[j]+fshift,(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 905 #else 906 PetscCall(PetscViewerASCIIPrintf(viewer,"%" PetscInt_FMT " %" PetscInt_FMT " %g\n", i+fshift, a->j[j]+fshift, (double)a->a[j])); 907 #endif 908 } 909 } 910 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 911 } else { 912 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_FALSE)); 913 if (A->factortype) { 914 for (i=0; i<m; i++) { 915 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 916 /* L part */ 917 for (j=a->i[i]; j<a->i[i+1]; j++) { 918 #if defined(PETSC_USE_COMPLEX) 919 if (PetscImaginaryPart(a->a[j]) > 0.0) { 920 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 921 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 922 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 923 } else { 924 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 925 } 926 #else 927 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 928 #endif 929 } 930 /* diagonal */ 931 j = a->diag[i]; 932 #if defined(PETSC_USE_COMPLEX) 933 if (PetscImaginaryPart(a->a[j]) > 0.0) { 934 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)PetscImaginaryPart(1.0/a->a[j]))); 935 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 936 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(1.0/a->a[j]),(double)(-PetscImaginaryPart(1.0/a->a[j])))); 937 } else { 938 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(1.0/a->a[j]))); 939 } 940 #else 941 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)(1.0/a->a[j]))); 942 #endif 943 944 /* U part */ 945 for (j=a->diag[i+1]+1; j<a->diag[i]; j++) { 946 #if defined(PETSC_USE_COMPLEX) 947 if (PetscImaginaryPart(a->a[j]) > 0.0) { 948 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 949 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 950 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)(-PetscImaginaryPart(a->a[j])))); 951 } else { 952 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 953 } 954 #else 955 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 956 #endif 957 } 958 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 959 } 960 } else { 961 for (i=0; i<m; i++) { 962 PetscCall(PetscViewerASCIIPrintf(viewer,"row %" PetscInt_FMT ":",i)); 963 for (j=a->i[i]; j<a->i[i+1]; j++) { 964 #if defined(PETSC_USE_COMPLEX) 965 if (PetscImaginaryPart(a->a[j]) > 0.0) { 966 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g + %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)PetscImaginaryPart(a->a[j]))); 967 } else if (PetscImaginaryPart(a->a[j]) < 0.0) { 968 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g - %g i)",a->j[j],(double)PetscRealPart(a->a[j]),(double)-PetscImaginaryPart(a->a[j]))); 969 } else { 970 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)PetscRealPart(a->a[j]))); 971 } 972 #else 973 PetscCall(PetscViewerASCIIPrintf(viewer," (%" PetscInt_FMT ", %g) ",a->j[j],(double)a->a[j])); 974 #endif 975 } 976 PetscCall(PetscViewerASCIIPrintf(viewer,"\n")); 977 } 978 } 979 PetscCall(PetscViewerASCIIUseTabs(viewer,PETSC_TRUE)); 980 } 981 PetscCall(PetscViewerFlush(viewer)); 982 PetscFunctionReturn(0); 983 } 984 985 #include <petscdraw.h> 986 PetscErrorCode MatView_SeqAIJ_Draw_Zoom(PetscDraw draw,void *Aa) 987 { 988 Mat A = (Mat) Aa; 989 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 990 PetscInt i,j,m = A->rmap->n; 991 int color; 992 PetscReal xl,yl,xr,yr,x_l,x_r,y_l,y_r; 993 PetscViewer viewer; 994 PetscViewerFormat format; 995 const PetscScalar *aa; 996 997 PetscFunctionBegin; 998 PetscCall(PetscObjectQuery((PetscObject)A,"Zoomviewer",(PetscObject*)&viewer)); 999 PetscCall(PetscViewerGetFormat(viewer,&format)); 1000 PetscCall(PetscDrawGetCoordinates(draw,&xl,&yl,&xr,&yr)); 1001 1002 /* loop over matrix elements drawing boxes */ 1003 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1004 if (format != PETSC_VIEWER_DRAW_CONTOUR) { 1005 PetscErrorCode ierr; 1006 ierr = PetscDrawCollectiveBegin(draw);PetscCall(ierr); 1007 /* Blue for negative, Cyan for zero and Red for positive */ 1008 color = PETSC_DRAW_BLUE; 1009 for (i=0; i<m; i++) { 1010 y_l = m - i - 1.0; y_r = y_l + 1.0; 1011 for (j=a->i[i]; j<a->i[i+1]; j++) { 1012 x_l = a->j[j]; x_r = x_l + 1.0; 1013 if (PetscRealPart(aa[j]) >= 0.) continue; 1014 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1015 } 1016 } 1017 color = PETSC_DRAW_CYAN; 1018 for (i=0; i<m; i++) { 1019 y_l = m - i - 1.0; y_r = y_l + 1.0; 1020 for (j=a->i[i]; j<a->i[i+1]; j++) { 1021 x_l = a->j[j]; x_r = x_l + 1.0; 1022 if (aa[j] != 0.) continue; 1023 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1024 } 1025 } 1026 color = PETSC_DRAW_RED; 1027 for (i=0; i<m; i++) { 1028 y_l = m - i - 1.0; y_r = y_l + 1.0; 1029 for (j=a->i[i]; j<a->i[i+1]; j++) { 1030 x_l = a->j[j]; x_r = x_l + 1.0; 1031 if (PetscRealPart(aa[j]) <= 0.) continue; 1032 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1033 } 1034 } 1035 ierr = PetscDrawCollectiveEnd(draw);PetscCall(ierr); 1036 } else { 1037 /* use contour shading to indicate magnitude of values */ 1038 /* first determine max of all nonzero values */ 1039 PetscReal minv = 0.0, maxv = 0.0; 1040 PetscInt nz = a->nz, count = 0; 1041 PetscDraw popup; 1042 PetscErrorCode ierr; 1043 1044 for (i=0; i<nz; i++) { 1045 if (PetscAbsScalar(aa[i]) > maxv) maxv = PetscAbsScalar(aa[i]); 1046 } 1047 if (minv >= maxv) maxv = minv + PETSC_SMALL; 1048 PetscCall(PetscDrawGetPopup(draw,&popup)); 1049 PetscCall(PetscDrawScalePopup(popup,minv,maxv)); 1050 1051 ierr = PetscDrawCollectiveBegin(draw);PetscCall(ierr); 1052 for (i=0; i<m; i++) { 1053 y_l = m - i - 1.0; 1054 y_r = y_l + 1.0; 1055 for (j=a->i[i]; j<a->i[i+1]; j++) { 1056 x_l = a->j[j]; 1057 x_r = x_l + 1.0; 1058 color = PetscDrawRealToColor(PetscAbsScalar(aa[count]),minv,maxv); 1059 PetscCall(PetscDrawRectangle(draw,x_l,y_l,x_r,y_r,color,color,color,color)); 1060 count++; 1061 } 1062 } 1063 ierr = PetscDrawCollectiveEnd(draw);PetscCall(ierr); 1064 } 1065 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1066 PetscFunctionReturn(0); 1067 } 1068 1069 #include <petscdraw.h> 1070 PetscErrorCode MatView_SeqAIJ_Draw(Mat A,PetscViewer viewer) 1071 { 1072 PetscDraw draw; 1073 PetscReal xr,yr,xl,yl,h,w; 1074 PetscBool isnull; 1075 1076 PetscFunctionBegin; 1077 PetscCall(PetscViewerDrawGetDraw(viewer,0,&draw)); 1078 PetscCall(PetscDrawIsNull(draw,&isnull)); 1079 if (isnull) PetscFunctionReturn(0); 1080 1081 xr = A->cmap->n; yr = A->rmap->n; h = yr/10.0; w = xr/10.0; 1082 xr += w; yr += h; xl = -w; yl = -h; 1083 PetscCall(PetscDrawSetCoordinates(draw,xl,yl,xr,yr)); 1084 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",(PetscObject)viewer)); 1085 PetscCall(PetscDrawZoom(draw,MatView_SeqAIJ_Draw_Zoom,A)); 1086 PetscCall(PetscObjectCompose((PetscObject)A,"Zoomviewer",NULL)); 1087 PetscCall(PetscDrawSave(draw)); 1088 PetscFunctionReturn(0); 1089 } 1090 1091 PetscErrorCode MatView_SeqAIJ(Mat A,PetscViewer viewer) 1092 { 1093 PetscBool iascii,isbinary,isdraw; 1094 1095 PetscFunctionBegin; 1096 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 1097 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 1098 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw)); 1099 if (iascii) { 1100 PetscCall(MatView_SeqAIJ_ASCII(A,viewer)); 1101 } else if (isbinary) { 1102 PetscCall(MatView_SeqAIJ_Binary(A,viewer)); 1103 } else if (isdraw) { 1104 PetscCall(MatView_SeqAIJ_Draw(A,viewer)); 1105 } 1106 PetscCall(MatView_SeqAIJ_Inode(A,viewer)); 1107 PetscFunctionReturn(0); 1108 } 1109 1110 PetscErrorCode MatAssemblyEnd_SeqAIJ(Mat A,MatAssemblyType mode) 1111 { 1112 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1113 PetscInt fshift = 0,i,*ai = a->i,*aj = a->j,*imax = a->imax; 1114 PetscInt m = A->rmap->n,*ip,N,*ailen = a->ilen,rmax = 0; 1115 MatScalar *aa = a->a,*ap; 1116 PetscReal ratio = 0.6; 1117 1118 PetscFunctionBegin; 1119 if (mode == MAT_FLUSH_ASSEMBLY) PetscFunctionReturn(0); 1120 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1121 if (A->was_assembled && A->ass_nonzerostate == A->nonzerostate) { 1122 /* we need to respect users asking to use or not the inodes routine in between matrix assemblies */ 1123 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1124 PetscFunctionReturn(0); 1125 } 1126 1127 if (m) rmax = ailen[0]; /* determine row with most nonzeros */ 1128 for (i=1; i<m; i++) { 1129 /* move each row back by the amount of empty slots (fshift) before it*/ 1130 fshift += imax[i-1] - ailen[i-1]; 1131 rmax = PetscMax(rmax,ailen[i]); 1132 if (fshift) { 1133 ip = aj + ai[i]; 1134 ap = aa + ai[i]; 1135 N = ailen[i]; 1136 PetscCall(PetscArraymove(ip-fshift,ip,N)); 1137 if (!A->structure_only) { 1138 PetscCall(PetscArraymove(ap-fshift,ap,N)); 1139 } 1140 } 1141 ai[i] = ai[i-1] + ailen[i-1]; 1142 } 1143 if (m) { 1144 fshift += imax[m-1] - ailen[m-1]; 1145 ai[m] = ai[m-1] + ailen[m-1]; 1146 } 1147 1148 /* reset ilen and imax for each row */ 1149 a->nonzerorowcnt = 0; 1150 if (A->structure_only) { 1151 PetscCall(PetscFree(a->imax)); 1152 PetscCall(PetscFree(a->ilen)); 1153 } else { /* !A->structure_only */ 1154 for (i=0; i<m; i++) { 1155 ailen[i] = imax[i] = ai[i+1] - ai[i]; 1156 a->nonzerorowcnt += ((ai[i+1] - ai[i]) > 0); 1157 } 1158 } 1159 a->nz = ai[m]; 1160 PetscCheckFalse(fshift && a->nounused == -1,PETSC_COMM_SELF,PETSC_ERR_PLIB, "Unused space detected in matrix: %" PetscInt_FMT " X %" PetscInt_FMT ", %" PetscInt_FMT " unneeded", m, A->cmap->n, fshift); 1161 1162 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1163 PetscCall(PetscInfo(A,"Matrix size: %" PetscInt_FMT " X %" PetscInt_FMT "; storage space: %" PetscInt_FMT " unneeded,%" PetscInt_FMT " used\n",m,A->cmap->n,fshift,a->nz)); 1164 PetscCall(PetscInfo(A,"Number of mallocs during MatSetValues() is %" PetscInt_FMT "\n",a->reallocs)); 1165 PetscCall(PetscInfo(A,"Maximum nonzeros in any row is %" PetscInt_FMT "\n",rmax)); 1166 1167 A->info.mallocs += a->reallocs; 1168 a->reallocs = 0; 1169 A->info.nz_unneeded = (PetscReal)fshift; 1170 a->rmax = rmax; 1171 1172 if (!A->structure_only) { 1173 PetscCall(MatCheckCompressedRow(A,a->nonzerorowcnt,&a->compressedrow,a->i,m,ratio)); 1174 } 1175 PetscCall(MatAssemblyEnd_SeqAIJ_Inode(A,mode)); 1176 PetscFunctionReturn(0); 1177 } 1178 1179 PetscErrorCode MatRealPart_SeqAIJ(Mat A) 1180 { 1181 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1182 PetscInt i,nz = a->nz; 1183 MatScalar *aa; 1184 1185 PetscFunctionBegin; 1186 PetscCall(MatSeqAIJGetArray(A,&aa)); 1187 for (i=0; i<nz; i++) aa[i] = PetscRealPart(aa[i]); 1188 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1189 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1190 PetscFunctionReturn(0); 1191 } 1192 1193 PetscErrorCode MatImaginaryPart_SeqAIJ(Mat A) 1194 { 1195 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1196 PetscInt i,nz = a->nz; 1197 MatScalar *aa; 1198 1199 PetscFunctionBegin; 1200 PetscCall(MatSeqAIJGetArray(A,&aa)); 1201 for (i=0; i<nz; i++) aa[i] = PetscImaginaryPart(aa[i]); 1202 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 1203 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1204 PetscFunctionReturn(0); 1205 } 1206 1207 PetscErrorCode MatZeroEntries_SeqAIJ(Mat A) 1208 { 1209 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1210 MatScalar *aa; 1211 1212 PetscFunctionBegin; 1213 PetscCall(MatSeqAIJGetArrayWrite(A,&aa)); 1214 PetscCall(PetscArrayzero(aa,a->i[A->rmap->n])); 1215 PetscCall(MatSeqAIJRestoreArrayWrite(A,&aa)); 1216 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 1217 PetscFunctionReturn(0); 1218 } 1219 1220 PETSC_INTERN PetscErrorCode MatResetPreallocationCOO_SeqAIJ(Mat A) 1221 { 1222 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1223 1224 PetscFunctionBegin; 1225 PetscCall(PetscFree(a->perm)); 1226 PetscCall(PetscFree(a->jmap)); 1227 PetscFunctionReturn(0); 1228 } 1229 1230 PetscErrorCode MatDestroy_SeqAIJ(Mat A) 1231 { 1232 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1233 1234 PetscFunctionBegin; 1235 #if defined(PETSC_USE_LOG) 1236 PetscLogObjectState((PetscObject)A,"Rows=%" PetscInt_FMT ", Cols=%" PetscInt_FMT ", NZ=%" PetscInt_FMT,A->rmap->n,A->cmap->n,a->nz); 1237 #endif 1238 PetscCall(MatSeqXAIJFreeAIJ(A,&a->a,&a->j,&a->i)); 1239 PetscCall(MatResetPreallocationCOO_SeqAIJ(A)); 1240 PetscCall(ISDestroy(&a->row)); 1241 PetscCall(ISDestroy(&a->col)); 1242 PetscCall(PetscFree(a->diag)); 1243 PetscCall(PetscFree(a->ibdiag)); 1244 PetscCall(PetscFree(a->imax)); 1245 PetscCall(PetscFree(a->ilen)); 1246 PetscCall(PetscFree(a->ipre)); 1247 PetscCall(PetscFree3(a->idiag,a->mdiag,a->ssor_work)); 1248 PetscCall(PetscFree(a->solve_work)); 1249 PetscCall(ISDestroy(&a->icol)); 1250 PetscCall(PetscFree(a->saved_values)); 1251 PetscCall(PetscFree2(a->compressedrow.i,a->compressedrow.rindex)); 1252 PetscCall(MatDestroy_SeqAIJ_Inode(A)); 1253 PetscCall(PetscFree(A->data)); 1254 1255 /* MatMatMultNumeric_SeqAIJ_SeqAIJ_Sorted may allocate this. 1256 That function is so heavily used (sometimes in an hidden way through multnumeric function pointers) 1257 that is hard to properly add this data to the MatProduct data. We free it here to avoid 1258 users reusing the matrix object with different data to incur in obscure segmentation faults 1259 due to different matrix sizes */ 1260 PetscCall(PetscObjectCompose((PetscObject)A,"__PETSc__ab_dense",NULL)); 1261 1262 PetscCall(PetscObjectChangeTypeName((PetscObject)A,NULL)); 1263 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetColumnIndices_C",NULL)); 1264 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatStoreValues_C",NULL)); 1265 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatRetrieveValues_C",NULL)); 1266 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsbaij_C",NULL)); 1267 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqbaij_C",NULL)); 1268 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijperm_C",NULL)); 1269 #if defined(PETSC_HAVE_CUDA) 1270 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcusparse_C",NULL)); 1271 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",NULL)); 1272 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",NULL)); 1273 #endif 1274 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 1275 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijkokkos_C",NULL)); 1276 #endif 1277 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqaijcrl_C",NULL)); 1278 #if defined(PETSC_HAVE_ELEMENTAL) 1279 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_elemental_C",NULL)); 1280 #endif 1281 #if defined(PETSC_HAVE_SCALAPACK) 1282 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_scalapack_C",NULL)); 1283 #endif 1284 #if defined(PETSC_HAVE_HYPRE) 1285 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_hypre_C",NULL)); 1286 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",NULL)); 1287 #endif 1288 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqdense_C",NULL)); 1289 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_seqsell_C",NULL)); 1290 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatConvert_seqaij_is_C",NULL)); 1291 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatIsTranspose_C",NULL)); 1292 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocation_C",NULL)); 1293 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatResetPreallocation_C",NULL)); 1294 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJSetPreallocationCSR_C",NULL)); 1295 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatReorderForNonzeroDiagonal_C",NULL)); 1296 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_is_seqaij_C",NULL)); 1297 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqdense_seqaij_C",NULL)); 1298 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatProductSetFromOptions_seqaij_seqaij_C",NULL)); 1299 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSeqAIJKron_C",NULL)); 1300 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetPreallocationCOO_C",NULL)); 1301 PetscCall(PetscObjectComposeFunction((PetscObject)A,"MatSetValuesCOO_C",NULL)); 1302 PetscFunctionReturn(0); 1303 } 1304 1305 PetscErrorCode MatSetOption_SeqAIJ(Mat A,MatOption op,PetscBool flg) 1306 { 1307 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1308 1309 PetscFunctionBegin; 1310 switch (op) { 1311 case MAT_ROW_ORIENTED: 1312 a->roworiented = flg; 1313 break; 1314 case MAT_KEEP_NONZERO_PATTERN: 1315 a->keepnonzeropattern = flg; 1316 break; 1317 case MAT_NEW_NONZERO_LOCATIONS: 1318 a->nonew = (flg ? 0 : 1); 1319 break; 1320 case MAT_NEW_NONZERO_LOCATION_ERR: 1321 a->nonew = (flg ? -1 : 0); 1322 break; 1323 case MAT_NEW_NONZERO_ALLOCATION_ERR: 1324 a->nonew = (flg ? -2 : 0); 1325 break; 1326 case MAT_UNUSED_NONZERO_LOCATION_ERR: 1327 a->nounused = (flg ? -1 : 0); 1328 break; 1329 case MAT_IGNORE_ZERO_ENTRIES: 1330 a->ignorezeroentries = flg; 1331 break; 1332 case MAT_SPD: 1333 case MAT_SYMMETRIC: 1334 case MAT_STRUCTURALLY_SYMMETRIC: 1335 case MAT_HERMITIAN: 1336 case MAT_SYMMETRY_ETERNAL: 1337 case MAT_STRUCTURE_ONLY: 1338 /* These options are handled directly by MatSetOption() */ 1339 break; 1340 case MAT_FORCE_DIAGONAL_ENTRIES: 1341 case MAT_IGNORE_OFF_PROC_ENTRIES: 1342 case MAT_USE_HASH_TABLE: 1343 PetscCall(PetscInfo(A,"Option %s ignored\n",MatOptions[op])); 1344 break; 1345 case MAT_USE_INODES: 1346 PetscCall(MatSetOption_SeqAIJ_Inode(A,MAT_USE_INODES,flg)); 1347 break; 1348 case MAT_SUBMAT_SINGLEIS: 1349 A->submat_singleis = flg; 1350 break; 1351 case MAT_SORTED_FULL: 1352 if (flg) A->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 1353 else A->ops->setvalues = MatSetValues_SeqAIJ; 1354 break; 1355 case MAT_FORM_EXPLICIT_TRANSPOSE: 1356 A->form_explicit_transpose = flg; 1357 break; 1358 default: 1359 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1360 } 1361 PetscFunctionReturn(0); 1362 } 1363 1364 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1365 { 1366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1367 PetscInt i,j,n,*ai=a->i,*aj=a->j; 1368 PetscScalar *x; 1369 const PetscScalar *aa; 1370 1371 PetscFunctionBegin; 1372 PetscCall(VecGetLocalSize(v,&n)); 1373 PetscCheckFalse(n != A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1374 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1375 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1376 PetscInt *diag=a->diag; 1377 PetscCall(VecGetArrayWrite(v,&x)); 1378 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1379 PetscCall(VecRestoreArrayWrite(v,&x)); 1380 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1381 PetscFunctionReturn(0); 1382 } 1383 1384 PetscCall(VecGetArrayWrite(v,&x)); 1385 for (i=0; i<n; i++) { 1386 x[i] = 0.0; 1387 for (j=ai[i]; j<ai[i+1]; j++) { 1388 if (aj[j] == i) { 1389 x[i] = aa[j]; 1390 break; 1391 } 1392 } 1393 } 1394 PetscCall(VecRestoreArrayWrite(v,&x)); 1395 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1396 PetscFunctionReturn(0); 1397 } 1398 1399 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1400 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1401 { 1402 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1403 const MatScalar *aa; 1404 PetscScalar *y; 1405 const PetscScalar *x; 1406 PetscInt m = A->rmap->n; 1407 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1408 const MatScalar *v; 1409 PetscScalar alpha; 1410 PetscInt n,i,j; 1411 const PetscInt *idx,*ii,*ridx=NULL; 1412 Mat_CompressedRow cprow = a->compressedrow; 1413 PetscBool usecprow = cprow.use; 1414 #endif 1415 1416 PetscFunctionBegin; 1417 if (zz != yy) PetscCall(VecCopy(zz,yy)); 1418 PetscCall(VecGetArrayRead(xx,&x)); 1419 PetscCall(VecGetArray(yy,&y)); 1420 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1421 1422 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1423 fortranmulttransposeaddaij_(&m,x,a->i,a->j,aa,y); 1424 #else 1425 if (usecprow) { 1426 m = cprow.nrows; 1427 ii = cprow.i; 1428 ridx = cprow.rindex; 1429 } else { 1430 ii = a->i; 1431 } 1432 for (i=0; i<m; i++) { 1433 idx = a->j + ii[i]; 1434 v = aa + ii[i]; 1435 n = ii[i+1] - ii[i]; 1436 if (usecprow) { 1437 alpha = x[ridx[i]]; 1438 } else { 1439 alpha = x[i]; 1440 } 1441 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1442 } 1443 #endif 1444 PetscCall(PetscLogFlops(2.0*a->nz)); 1445 PetscCall(VecRestoreArrayRead(xx,&x)); 1446 PetscCall(VecRestoreArray(yy,&y)); 1447 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1448 PetscFunctionReturn(0); 1449 } 1450 1451 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1452 { 1453 PetscFunctionBegin; 1454 PetscCall(VecSet(yy,0.0)); 1455 PetscCall(MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy)); 1456 PetscFunctionReturn(0); 1457 } 1458 1459 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1460 1461 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1462 { 1463 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1464 PetscScalar *y; 1465 const PetscScalar *x; 1466 const MatScalar *aa,*a_a; 1467 PetscInt m=A->rmap->n; 1468 const PetscInt *aj,*ii,*ridx=NULL; 1469 PetscInt n,i; 1470 PetscScalar sum; 1471 PetscBool usecprow=a->compressedrow.use; 1472 1473 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1474 #pragma disjoint(*x,*y,*aa) 1475 #endif 1476 1477 PetscFunctionBegin; 1478 if (a->inode.use && a->inode.checked) { 1479 PetscCall(MatMult_SeqAIJ_Inode(A,xx,yy)); 1480 PetscFunctionReturn(0); 1481 } 1482 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1483 PetscCall(VecGetArrayRead(xx,&x)); 1484 PetscCall(VecGetArray(yy,&y)); 1485 ii = a->i; 1486 if (usecprow) { /* use compressed row format */ 1487 PetscCall(PetscArrayzero(y,m)); 1488 m = a->compressedrow.nrows; 1489 ii = a->compressedrow.i; 1490 ridx = a->compressedrow.rindex; 1491 for (i=0; i<m; i++) { 1492 n = ii[i+1] - ii[i]; 1493 aj = a->j + ii[i]; 1494 aa = a_a + ii[i]; 1495 sum = 0.0; 1496 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1497 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1498 y[*ridx++] = sum; 1499 } 1500 } else { /* do not use compressed row format */ 1501 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1502 aj = a->j; 1503 aa = a_a; 1504 fortranmultaij_(&m,x,ii,aj,aa,y); 1505 #else 1506 for (i=0; i<m; i++) { 1507 n = ii[i+1] - ii[i]; 1508 aj = a->j + ii[i]; 1509 aa = a_a + ii[i]; 1510 sum = 0.0; 1511 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1512 y[i] = sum; 1513 } 1514 #endif 1515 } 1516 PetscCall(PetscLogFlops(2.0*a->nz - a->nonzerorowcnt)); 1517 PetscCall(VecRestoreArrayRead(xx,&x)); 1518 PetscCall(VecRestoreArray(yy,&y)); 1519 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1520 PetscFunctionReturn(0); 1521 } 1522 1523 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1524 { 1525 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1526 PetscScalar *y; 1527 const PetscScalar *x; 1528 const MatScalar *aa,*a_a; 1529 PetscInt m=A->rmap->n; 1530 const PetscInt *aj,*ii,*ridx=NULL; 1531 PetscInt n,i,nonzerorow=0; 1532 PetscScalar sum; 1533 PetscBool usecprow=a->compressedrow.use; 1534 1535 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1536 #pragma disjoint(*x,*y,*aa) 1537 #endif 1538 1539 PetscFunctionBegin; 1540 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1541 PetscCall(VecGetArrayRead(xx,&x)); 1542 PetscCall(VecGetArray(yy,&y)); 1543 if (usecprow) { /* use compressed row format */ 1544 m = a->compressedrow.nrows; 1545 ii = a->compressedrow.i; 1546 ridx = a->compressedrow.rindex; 1547 for (i=0; i<m; i++) { 1548 n = ii[i+1] - ii[i]; 1549 aj = a->j + ii[i]; 1550 aa = a_a + ii[i]; 1551 sum = 0.0; 1552 nonzerorow += (n>0); 1553 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1554 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1555 y[*ridx++] = sum; 1556 } 1557 } else { /* do not use compressed row format */ 1558 ii = a->i; 1559 for (i=0; i<m; i++) { 1560 n = ii[i+1] - ii[i]; 1561 aj = a->j + ii[i]; 1562 aa = a_a + ii[i]; 1563 sum = 0.0; 1564 nonzerorow += (n>0); 1565 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1566 y[i] = sum; 1567 } 1568 } 1569 PetscCall(PetscLogFlops(2.0*a->nz - nonzerorow)); 1570 PetscCall(VecRestoreArrayRead(xx,&x)); 1571 PetscCall(VecRestoreArray(yy,&y)); 1572 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1573 PetscFunctionReturn(0); 1574 } 1575 1576 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1577 { 1578 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1579 PetscScalar *y,*z; 1580 const PetscScalar *x; 1581 const MatScalar *aa,*a_a; 1582 PetscInt m = A->rmap->n,*aj,*ii; 1583 PetscInt n,i,*ridx=NULL; 1584 PetscScalar sum; 1585 PetscBool usecprow=a->compressedrow.use; 1586 1587 PetscFunctionBegin; 1588 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1589 PetscCall(VecGetArrayRead(xx,&x)); 1590 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1591 if (usecprow) { /* use compressed row format */ 1592 if (zz != yy) { 1593 PetscCall(PetscArraycpy(z,y,m)); 1594 } 1595 m = a->compressedrow.nrows; 1596 ii = a->compressedrow.i; 1597 ridx = a->compressedrow.rindex; 1598 for (i=0; i<m; i++) { 1599 n = ii[i+1] - ii[i]; 1600 aj = a->j + ii[i]; 1601 aa = a_a + ii[i]; 1602 sum = y[*ridx]; 1603 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1604 z[*ridx++] = sum; 1605 } 1606 } else { /* do not use compressed row format */ 1607 ii = a->i; 1608 for (i=0; i<m; i++) { 1609 n = ii[i+1] - ii[i]; 1610 aj = a->j + ii[i]; 1611 aa = a_a + ii[i]; 1612 sum = y[i]; 1613 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1614 z[i] = sum; 1615 } 1616 } 1617 PetscCall(PetscLogFlops(2.0*a->nz)); 1618 PetscCall(VecRestoreArrayRead(xx,&x)); 1619 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1620 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1621 PetscFunctionReturn(0); 1622 } 1623 1624 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1625 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1626 { 1627 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1628 PetscScalar *y,*z; 1629 const PetscScalar *x; 1630 const MatScalar *aa,*a_a; 1631 const PetscInt *aj,*ii,*ridx=NULL; 1632 PetscInt m = A->rmap->n,n,i; 1633 PetscScalar sum; 1634 PetscBool usecprow=a->compressedrow.use; 1635 1636 PetscFunctionBegin; 1637 if (a->inode.use && a->inode.checked) { 1638 PetscCall(MatMultAdd_SeqAIJ_Inode(A,xx,yy,zz)); 1639 PetscFunctionReturn(0); 1640 } 1641 PetscCall(MatSeqAIJGetArrayRead(A,&a_a)); 1642 PetscCall(VecGetArrayRead(xx,&x)); 1643 PetscCall(VecGetArrayPair(yy,zz,&y,&z)); 1644 if (usecprow) { /* use compressed row format */ 1645 if (zz != yy) { 1646 PetscCall(PetscArraycpy(z,y,m)); 1647 } 1648 m = a->compressedrow.nrows; 1649 ii = a->compressedrow.i; 1650 ridx = a->compressedrow.rindex; 1651 for (i=0; i<m; i++) { 1652 n = ii[i+1] - ii[i]; 1653 aj = a->j + ii[i]; 1654 aa = a_a + ii[i]; 1655 sum = y[*ridx]; 1656 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1657 z[*ridx++] = sum; 1658 } 1659 } else { /* do not use compressed row format */ 1660 ii = a->i; 1661 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1662 aj = a->j; 1663 aa = a_a; 1664 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1665 #else 1666 for (i=0; i<m; i++) { 1667 n = ii[i+1] - ii[i]; 1668 aj = a->j + ii[i]; 1669 aa = a_a + ii[i]; 1670 sum = y[i]; 1671 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1672 z[i] = sum; 1673 } 1674 #endif 1675 } 1676 PetscCall(PetscLogFlops(2.0*a->nz)); 1677 PetscCall(VecRestoreArrayRead(xx,&x)); 1678 PetscCall(VecRestoreArrayPair(yy,zz,&y,&z)); 1679 PetscCall(MatSeqAIJRestoreArrayRead(A,&a_a)); 1680 PetscFunctionReturn(0); 1681 } 1682 1683 /* 1684 Adds diagonal pointers to sparse matrix structure. 1685 */ 1686 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1687 { 1688 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1689 PetscInt i,j,m = A->rmap->n; 1690 PetscBool alreadySet = PETSC_TRUE; 1691 1692 PetscFunctionBegin; 1693 if (!a->diag) { 1694 PetscCall(PetscMalloc1(m,&a->diag)); 1695 PetscCall(PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt))); 1696 alreadySet = PETSC_FALSE; 1697 } 1698 for (i=0; i<A->rmap->n; i++) { 1699 /* If A's diagonal is already correctly set, this fast track enables cheap and repeated MatMarkDiagonal_SeqAIJ() calls */ 1700 if (alreadySet) { 1701 PetscInt pos = a->diag[i]; 1702 if (pos >= a->i[i] && pos < a->i[i+1] && a->j[pos] == i) continue; 1703 } 1704 1705 a->diag[i] = a->i[i+1]; 1706 for (j=a->i[i]; j<a->i[i+1]; j++) { 1707 if (a->j[j] == i) { 1708 a->diag[i] = j; 1709 break; 1710 } 1711 } 1712 } 1713 PetscFunctionReturn(0); 1714 } 1715 1716 PetscErrorCode MatShift_SeqAIJ(Mat A,PetscScalar v) 1717 { 1718 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1719 const PetscInt *diag = (const PetscInt*)a->diag; 1720 const PetscInt *ii = (const PetscInt*) a->i; 1721 PetscInt i,*mdiag = NULL; 1722 PetscInt cnt = 0; /* how many diagonals are missing */ 1723 1724 PetscFunctionBegin; 1725 if (!A->preallocated || !a->nz) { 1726 PetscCall(MatSeqAIJSetPreallocation(A,1,NULL)); 1727 PetscCall(MatShift_Basic(A,v)); 1728 PetscFunctionReturn(0); 1729 } 1730 1731 if (a->diagonaldense) { 1732 cnt = 0; 1733 } else { 1734 PetscCall(PetscCalloc1(A->rmap->n,&mdiag)); 1735 for (i=0; i<A->rmap->n; i++) { 1736 if (i < A->cmap->n && diag[i] >= ii[i+1]) { /* 'out of range' rows never have diagonals */ 1737 cnt++; 1738 mdiag[i] = 1; 1739 } 1740 } 1741 } 1742 if (!cnt) { 1743 PetscCall(MatShift_Basic(A,v)); 1744 } else { 1745 PetscScalar *olda = a->a; /* preserve pointers to current matrix nonzeros structure and values */ 1746 PetscInt *oldj = a->j, *oldi = a->i; 1747 PetscBool singlemalloc = a->singlemalloc,free_a = a->free_a,free_ij = a->free_ij; 1748 1749 a->a = NULL; 1750 a->j = NULL; 1751 a->i = NULL; 1752 /* increase the values in imax for each row where a diagonal is being inserted then reallocate the matrix data structures */ 1753 for (i=0; i<PetscMin(A->rmap->n,A->cmap->n); i++) { 1754 a->imax[i] += mdiag[i]; 1755 } 1756 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,0,a->imax)); 1757 1758 /* copy old values into new matrix data structure */ 1759 for (i=0; i<A->rmap->n; i++) { 1760 PetscCall(MatSetValues(A,1,&i,a->imax[i] - mdiag[i],&oldj[oldi[i]],&olda[oldi[i]],ADD_VALUES)); 1761 if (i < A->cmap->n) { 1762 PetscCall(MatSetValue(A,i,i,v,ADD_VALUES)); 1763 } 1764 } 1765 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1766 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 1767 if (singlemalloc) { 1768 PetscCall(PetscFree3(olda,oldj,oldi)); 1769 } else { 1770 if (free_a) PetscCall(PetscFree(olda)); 1771 if (free_ij) PetscCall(PetscFree(oldj)); 1772 if (free_ij) PetscCall(PetscFree(oldi)); 1773 } 1774 } 1775 PetscCall(PetscFree(mdiag)); 1776 a->diagonaldense = PETSC_TRUE; 1777 PetscFunctionReturn(0); 1778 } 1779 1780 /* 1781 Checks for missing diagonals 1782 */ 1783 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1784 { 1785 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1786 PetscInt *diag,*ii = a->i,i; 1787 1788 PetscFunctionBegin; 1789 *missing = PETSC_FALSE; 1790 if (A->rmap->n > 0 && !ii) { 1791 *missing = PETSC_TRUE; 1792 if (d) *d = 0; 1793 PetscCall(PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n")); 1794 } else { 1795 PetscInt n; 1796 n = PetscMin(A->rmap->n, A->cmap->n); 1797 diag = a->diag; 1798 for (i=0; i<n; i++) { 1799 if (diag[i] >= ii[i+1]) { 1800 *missing = PETSC_TRUE; 1801 if (d) *d = i; 1802 PetscCall(PetscInfo(A,"Matrix is missing diagonal number %" PetscInt_FMT "\n",i)); 1803 break; 1804 } 1805 } 1806 } 1807 PetscFunctionReturn(0); 1808 } 1809 1810 #include <petscblaslapack.h> 1811 #include <petsc/private/kernels/blockinvert.h> 1812 1813 /* 1814 Note that values is allocated externally by the PC and then passed into this routine 1815 */ 1816 PetscErrorCode MatInvertVariableBlockDiagonal_SeqAIJ(Mat A,PetscInt nblocks,const PetscInt *bsizes,PetscScalar *diag) 1817 { 1818 PetscInt n = A->rmap->n, i, ncnt = 0, *indx,j,bsizemax = 0,*v_pivots; 1819 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 1820 const PetscReal shift = 0.0; 1821 PetscInt ipvt[5]; 1822 PetscScalar work[25],*v_work; 1823 1824 PetscFunctionBegin; 1825 allowzeropivot = PetscNot(A->erroriffailure); 1826 for (i=0; i<nblocks; i++) ncnt += bsizes[i]; 1827 PetscCheckFalse(ncnt != n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Total blocksizes %" PetscInt_FMT " doesn't match number matrix rows %" PetscInt_FMT,ncnt,n); 1828 for (i=0; i<nblocks; i++) { 1829 bsizemax = PetscMax(bsizemax,bsizes[i]); 1830 } 1831 PetscCall(PetscMalloc1(bsizemax,&indx)); 1832 if (bsizemax > 7) { 1833 PetscCall(PetscMalloc2(bsizemax,&v_work,bsizemax,&v_pivots)); 1834 } 1835 ncnt = 0; 1836 for (i=0; i<nblocks; i++) { 1837 for (j=0; j<bsizes[i]; j++) indx[j] = ncnt+j; 1838 PetscCall(MatGetValues(A,bsizes[i],indx,bsizes[i],indx,diag)); 1839 switch (bsizes[i]) { 1840 case 1: 1841 *diag = 1.0/(*diag); 1842 break; 1843 case 2: 1844 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 1845 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1846 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 1847 break; 1848 case 3: 1849 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 1850 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1851 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 1852 break; 1853 case 4: 1854 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 1855 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1856 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 1857 break; 1858 case 5: 1859 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 1860 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1861 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 1862 break; 1863 case 6: 1864 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 1865 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1866 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 1867 break; 1868 case 7: 1869 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 1870 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1871 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 1872 break; 1873 default: 1874 PetscCall(PetscKernel_A_gets_inverse_A(bsizes[i],diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 1875 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1876 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bsizes[i])); 1877 } 1878 ncnt += bsizes[i]; 1879 diag += bsizes[i]*bsizes[i]; 1880 } 1881 if (bsizemax > 7) { 1882 PetscCall(PetscFree2(v_work,v_pivots)); 1883 } 1884 PetscCall(PetscFree(indx)); 1885 PetscFunctionReturn(0); 1886 } 1887 1888 /* 1889 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1890 */ 1891 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1892 { 1893 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1894 PetscInt i,*diag,m = A->rmap->n; 1895 const MatScalar *v; 1896 PetscScalar *idiag,*mdiag; 1897 1898 PetscFunctionBegin; 1899 if (a->idiagvalid) PetscFunctionReturn(0); 1900 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 1901 diag = a->diag; 1902 if (!a->idiag) { 1903 PetscCall(PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work)); 1904 PetscCall(PetscLogObjectMemory((PetscObject)A,3*m*sizeof(PetscScalar))); 1905 } 1906 1907 mdiag = a->mdiag; 1908 idiag = a->idiag; 1909 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 1910 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1911 for (i=0; i<m; i++) { 1912 mdiag[i] = v[diag[i]]; 1913 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1914 if (PetscRealPart(fshift)) { 1915 PetscCall(PetscInfo(A,"Zero diagonal on row %" PetscInt_FMT "\n",i)); 1916 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1917 A->factorerror_zeropivot_value = 0.0; 1918 A->factorerror_zeropivot_row = i; 1919 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %" PetscInt_FMT,i); 1920 } 1921 idiag[i] = 1.0/v[diag[i]]; 1922 } 1923 PetscCall(PetscLogFlops(m)); 1924 } else { 1925 for (i=0; i<m; i++) { 1926 mdiag[i] = v[diag[i]]; 1927 idiag[i] = omega/(fshift + v[diag[i]]); 1928 } 1929 PetscCall(PetscLogFlops(2.0*m)); 1930 } 1931 a->idiagvalid = PETSC_TRUE; 1932 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 1933 PetscFunctionReturn(0); 1934 } 1935 1936 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1937 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1938 { 1939 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1940 PetscScalar *x,d,sum,*t,scale; 1941 const MatScalar *v,*idiag=NULL,*mdiag,*aa; 1942 const PetscScalar *b, *bs,*xb, *ts; 1943 PetscInt n,m = A->rmap->n,i; 1944 const PetscInt *idx,*diag; 1945 1946 PetscFunctionBegin; 1947 if (a->inode.use && a->inode.checked && omega == 1.0 && fshift == 0.0) { 1948 PetscCall(MatSOR_SeqAIJ_Inode(A,bb,omega,flag,fshift,its,lits,xx)); 1949 PetscFunctionReturn(0); 1950 } 1951 its = its*lits; 1952 1953 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1954 if (!a->idiagvalid) PetscCall(MatInvertDiagonal_SeqAIJ(A,omega,fshift)); 1955 a->fshift = fshift; 1956 a->omega = omega; 1957 1958 diag = a->diag; 1959 t = a->ssor_work; 1960 idiag = a->idiag; 1961 mdiag = a->mdiag; 1962 1963 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 1964 PetscCall(VecGetArray(xx,&x)); 1965 PetscCall(VecGetArrayRead(bb,&b)); 1966 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1967 if (flag == SOR_APPLY_UPPER) { 1968 /* apply (U + D/omega) to the vector */ 1969 bs = b; 1970 for (i=0; i<m; i++) { 1971 d = fshift + mdiag[i]; 1972 n = a->i[i+1] - diag[i] - 1; 1973 idx = a->j + diag[i] + 1; 1974 v = aa + diag[i] + 1; 1975 sum = b[i]*d/omega; 1976 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1977 x[i] = sum; 1978 } 1979 PetscCall(VecRestoreArray(xx,&x)); 1980 PetscCall(VecRestoreArrayRead(bb,&b)); 1981 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 1982 PetscCall(PetscLogFlops(a->nz)); 1983 PetscFunctionReturn(0); 1984 } 1985 1986 PetscCheckFalse(flag == SOR_APPLY_LOWER,PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1987 else if (flag & SOR_EISENSTAT) { 1988 /* Let A = L + U + D; where L is lower triangular, 1989 U is upper triangular, E = D/omega; This routine applies 1990 1991 (L + E)^{-1} A (U + E)^{-1} 1992 1993 to a vector efficiently using Eisenstat's trick. 1994 */ 1995 scale = (2.0/omega) - 1.0; 1996 1997 /* x = (E + U)^{-1} b */ 1998 for (i=m-1; i>=0; i--) { 1999 n = a->i[i+1] - diag[i] - 1; 2000 idx = a->j + diag[i] + 1; 2001 v = aa + diag[i] + 1; 2002 sum = b[i]; 2003 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2004 x[i] = sum*idiag[i]; 2005 } 2006 2007 /* t = b - (2*E - D)x */ 2008 v = aa; 2009 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 2010 2011 /* t = (E + L)^{-1}t */ 2012 ts = t; 2013 diag = a->diag; 2014 for (i=0; i<m; i++) { 2015 n = diag[i] - a->i[i]; 2016 idx = a->j + a->i[i]; 2017 v = aa + a->i[i]; 2018 sum = t[i]; 2019 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 2020 t[i] = sum*idiag[i]; 2021 /* x = x + t */ 2022 x[i] += t[i]; 2023 } 2024 2025 PetscCall(PetscLogFlops(6.0*m-1 + 2.0*a->nz)); 2026 PetscCall(VecRestoreArray(xx,&x)); 2027 PetscCall(VecRestoreArrayRead(bb,&b)); 2028 PetscFunctionReturn(0); 2029 } 2030 if (flag & SOR_ZERO_INITIAL_GUESS) { 2031 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2032 for (i=0; i<m; i++) { 2033 n = diag[i] - a->i[i]; 2034 idx = a->j + a->i[i]; 2035 v = aa + a->i[i]; 2036 sum = b[i]; 2037 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2038 t[i] = sum; 2039 x[i] = sum*idiag[i]; 2040 } 2041 xb = t; 2042 PetscCall(PetscLogFlops(a->nz)); 2043 } else xb = b; 2044 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2045 for (i=m-1; i>=0; i--) { 2046 n = a->i[i+1] - diag[i] - 1; 2047 idx = a->j + diag[i] + 1; 2048 v = aa + diag[i] + 1; 2049 sum = xb[i]; 2050 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2051 if (xb == b) { 2052 x[i] = sum*idiag[i]; 2053 } else { 2054 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2055 } 2056 } 2057 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2058 } 2059 its--; 2060 } 2061 while (its--) { 2062 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 2063 for (i=0; i<m; i++) { 2064 /* lower */ 2065 n = diag[i] - a->i[i]; 2066 idx = a->j + a->i[i]; 2067 v = aa + a->i[i]; 2068 sum = b[i]; 2069 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2070 t[i] = sum; /* save application of the lower-triangular part */ 2071 /* upper */ 2072 n = a->i[i+1] - diag[i] - 1; 2073 idx = a->j + diag[i] + 1; 2074 v = aa + diag[i] + 1; 2075 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2076 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2077 } 2078 xb = t; 2079 PetscCall(PetscLogFlops(2.0*a->nz)); 2080 } else xb = b; 2081 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 2082 for (i=m-1; i>=0; i--) { 2083 sum = xb[i]; 2084 if (xb == b) { 2085 /* whole matrix (no checkpointing available) */ 2086 n = a->i[i+1] - a->i[i]; 2087 idx = a->j + a->i[i]; 2088 v = aa + a->i[i]; 2089 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2090 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 2091 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 2092 n = a->i[i+1] - diag[i] - 1; 2093 idx = a->j + diag[i] + 1; 2094 v = aa + diag[i] + 1; 2095 PetscSparseDenseMinusDot(sum,x,v,idx,n); 2096 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 2097 } 2098 } 2099 if (xb == b) { 2100 PetscCall(PetscLogFlops(2.0*a->nz)); 2101 } else { 2102 PetscCall(PetscLogFlops(a->nz)); /* assumes 1/2 in upper */ 2103 } 2104 } 2105 } 2106 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2107 PetscCall(VecRestoreArray(xx,&x)); 2108 PetscCall(VecRestoreArrayRead(bb,&b)); 2109 PetscFunctionReturn(0); 2110 } 2111 2112 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 2113 { 2114 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2115 2116 PetscFunctionBegin; 2117 info->block_size = 1.0; 2118 info->nz_allocated = a->maxnz; 2119 info->nz_used = a->nz; 2120 info->nz_unneeded = (a->maxnz - a->nz); 2121 info->assemblies = A->num_ass; 2122 info->mallocs = A->info.mallocs; 2123 info->memory = ((PetscObject)A)->mem; 2124 if (A->factortype) { 2125 info->fill_ratio_given = A->info.fill_ratio_given; 2126 info->fill_ratio_needed = A->info.fill_ratio_needed; 2127 info->factor_mallocs = A->info.factor_mallocs; 2128 } else { 2129 info->fill_ratio_given = 0; 2130 info->fill_ratio_needed = 0; 2131 info->factor_mallocs = 0; 2132 } 2133 PetscFunctionReturn(0); 2134 } 2135 2136 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2137 { 2138 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2139 PetscInt i,m = A->rmap->n - 1; 2140 const PetscScalar *xx; 2141 PetscScalar *bb,*aa; 2142 PetscInt d = 0; 2143 2144 PetscFunctionBegin; 2145 if (x && b) { 2146 PetscCall(VecGetArrayRead(x,&xx)); 2147 PetscCall(VecGetArray(b,&bb)); 2148 for (i=0; i<N; i++) { 2149 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2150 if (rows[i] >= A->cmap->n) continue; 2151 bb[rows[i]] = diag*xx[rows[i]]; 2152 } 2153 PetscCall(VecRestoreArrayRead(x,&xx)); 2154 PetscCall(VecRestoreArray(b,&bb)); 2155 } 2156 2157 PetscCall(MatSeqAIJGetArray(A,&aa)); 2158 if (a->keepnonzeropattern) { 2159 for (i=0; i<N; i++) { 2160 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2161 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2162 } 2163 if (diag != 0.0) { 2164 for (i=0; i<N; i++) { 2165 d = rows[i]; 2166 if (rows[i] >= A->cmap->n) continue; 2167 PetscCheckFalse(a->diag[d] >= a->i[d+1],PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in the zeroed row %" PetscInt_FMT,d); 2168 } 2169 for (i=0; i<N; i++) { 2170 if (rows[i] >= A->cmap->n) continue; 2171 aa[a->diag[rows[i]]] = diag; 2172 } 2173 } 2174 } else { 2175 if (diag != 0.0) { 2176 for (i=0; i<N; i++) { 2177 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2178 if (a->ilen[rows[i]] > 0) { 2179 if (rows[i] >= A->cmap->n) { 2180 a->ilen[rows[i]] = 0; 2181 } else { 2182 a->ilen[rows[i]] = 1; 2183 aa[a->i[rows[i]]] = diag; 2184 a->j[a->i[rows[i]]] = rows[i]; 2185 } 2186 } else if (rows[i] < A->cmap->n) { /* in case row was completely empty */ 2187 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2188 } 2189 } 2190 } else { 2191 for (i=0; i<N; i++) { 2192 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2193 a->ilen[rows[i]] = 0; 2194 } 2195 } 2196 A->nonzerostate++; 2197 } 2198 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2199 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2200 PetscFunctionReturn(0); 2201 } 2202 2203 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 2204 { 2205 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2206 PetscInt i,j,m = A->rmap->n - 1,d = 0; 2207 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 2208 const PetscScalar *xx; 2209 PetscScalar *bb,*aa; 2210 2211 PetscFunctionBegin; 2212 if (!N) PetscFunctionReturn(0); 2213 PetscCall(MatSeqAIJGetArray(A,&aa)); 2214 if (x && b) { 2215 PetscCall(VecGetArrayRead(x,&xx)); 2216 PetscCall(VecGetArray(b,&bb)); 2217 vecs = PETSC_TRUE; 2218 } 2219 PetscCall(PetscCalloc1(A->rmap->n,&zeroed)); 2220 for (i=0; i<N; i++) { 2221 PetscCheckFalse(rows[i] < 0 || rows[i] > m,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %" PetscInt_FMT " out of range", rows[i]); 2222 PetscCall(PetscArrayzero(&aa[a->i[rows[i]]],a->ilen[rows[i]])); 2223 2224 zeroed[rows[i]] = PETSC_TRUE; 2225 } 2226 for (i=0; i<A->rmap->n; i++) { 2227 if (!zeroed[i]) { 2228 for (j=a->i[i]; j<a->i[i+1]; j++) { 2229 if (a->j[j] < A->rmap->n && zeroed[a->j[j]]) { 2230 if (vecs) bb[i] -= aa[j]*xx[a->j[j]]; 2231 aa[j] = 0.0; 2232 } 2233 } 2234 } else if (vecs && i < A->cmap->N) bb[i] = diag*xx[i]; 2235 } 2236 if (x && b) { 2237 PetscCall(VecRestoreArrayRead(x,&xx)); 2238 PetscCall(VecRestoreArray(b,&bb)); 2239 } 2240 PetscCall(PetscFree(zeroed)); 2241 if (diag != 0.0) { 2242 PetscCall(MatMissingDiagonal_SeqAIJ(A,&missing,&d)); 2243 if (missing) { 2244 for (i=0; i<N; i++) { 2245 if (rows[i] >= A->cmap->N) continue; 2246 PetscCheckFalse(a->nonew && rows[i] >= d,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")",d,rows[i]); 2247 PetscCall(MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES)); 2248 } 2249 } else { 2250 for (i=0; i<N; i++) { 2251 aa[a->diag[rows[i]]] = diag; 2252 } 2253 } 2254 } 2255 PetscCall(MatSeqAIJRestoreArray(A,&aa)); 2256 PetscCall((*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY)); 2257 PetscFunctionReturn(0); 2258 } 2259 2260 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2261 { 2262 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2263 const PetscScalar *aa; 2264 PetscInt *itmp; 2265 2266 PetscFunctionBegin; 2267 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2268 *nz = a->i[row+1] - a->i[row]; 2269 if (v) *v = (PetscScalar*)(aa + a->i[row]); 2270 if (idx) { 2271 itmp = a->j + a->i[row]; 2272 if (*nz) *idx = itmp; 2273 else *idx = NULL; 2274 } 2275 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2276 PetscFunctionReturn(0); 2277 } 2278 2279 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 2280 { 2281 PetscFunctionBegin; 2282 if (nz) *nz = 0; 2283 if (idx) *idx = NULL; 2284 if (v) *v = NULL; 2285 PetscFunctionReturn(0); 2286 } 2287 2288 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 2289 { 2290 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2291 const MatScalar *v; 2292 PetscReal sum = 0.0; 2293 PetscInt i,j; 2294 2295 PetscFunctionBegin; 2296 PetscCall(MatSeqAIJGetArrayRead(A,&v)); 2297 if (type == NORM_FROBENIUS) { 2298 #if defined(PETSC_USE_REAL___FP16) 2299 PetscBLASInt one = 1,nz = a->nz; 2300 PetscStackCallBLAS("BLASnrm2",*nrm = BLASnrm2_(&nz,v,&one)); 2301 #else 2302 for (i=0; i<a->nz; i++) { 2303 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 2304 } 2305 *nrm = PetscSqrtReal(sum); 2306 #endif 2307 PetscCall(PetscLogFlops(2.0*a->nz)); 2308 } else if (type == NORM_1) { 2309 PetscReal *tmp; 2310 PetscInt *jj = a->j; 2311 PetscCall(PetscCalloc1(A->cmap->n+1,&tmp)); 2312 *nrm = 0.0; 2313 for (j=0; j<a->nz; j++) { 2314 tmp[*jj++] += PetscAbsScalar(*v); v++; 2315 } 2316 for (j=0; j<A->cmap->n; j++) { 2317 if (tmp[j] > *nrm) *nrm = tmp[j]; 2318 } 2319 PetscCall(PetscFree(tmp)); 2320 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2321 } else if (type == NORM_INFINITY) { 2322 *nrm = 0.0; 2323 for (j=0; j<A->rmap->n; j++) { 2324 const PetscScalar *v2 = v + a->i[j]; 2325 sum = 0.0; 2326 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 2327 sum += PetscAbsScalar(*v2); v2++; 2328 } 2329 if (sum > *nrm) *nrm = sum; 2330 } 2331 PetscCall(PetscLogFlops(PetscMax(a->nz-1,0))); 2332 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 2333 PetscCall(MatSeqAIJRestoreArrayRead(A,&v)); 2334 PetscFunctionReturn(0); 2335 } 2336 2337 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 2338 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 2339 { 2340 PetscInt i,j,anzj; 2341 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 2342 PetscInt an=A->cmap->N,am=A->rmap->N; 2343 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 2344 2345 PetscFunctionBegin; 2346 /* Allocate space for symbolic transpose info and work array */ 2347 PetscCall(PetscCalloc1(an+1,&ati)); 2348 PetscCall(PetscMalloc1(ai[am],&atj)); 2349 PetscCall(PetscMalloc1(an,&atfill)); 2350 2351 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 2352 /* Note: offset by 1 for fast conversion into csr format. */ 2353 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2354 /* Form ati for csr format of A^T. */ 2355 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2356 2357 /* Copy ati into atfill so we have locations of the next free space in atj */ 2358 PetscCall(PetscArraycpy(atfill,ati,an)); 2359 2360 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2361 for (i=0;i<am;i++) { 2362 anzj = ai[i+1] - ai[i]; 2363 for (j=0;j<anzj;j++) { 2364 atj[atfill[*aj]] = i; 2365 atfill[*aj++] += 1; 2366 } 2367 } 2368 2369 /* Clean up temporary space and complete requests. */ 2370 PetscCall(PetscFree(atfill)); 2371 PetscCall(MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B)); 2372 PetscCall(MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs))); 2373 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2374 2375 b = (Mat_SeqAIJ*)((*B)->data); 2376 b->free_a = PETSC_FALSE; 2377 b->free_ij = PETSC_TRUE; 2378 b->nonew = 0; 2379 PetscFunctionReturn(0); 2380 } 2381 2382 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2383 { 2384 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2385 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2386 const MatScalar *va,*vb; 2387 PetscInt ma,na,mb,nb, i; 2388 2389 PetscFunctionBegin; 2390 PetscCall(MatGetSize(A,&ma,&na)); 2391 PetscCall(MatGetSize(B,&mb,&nb)); 2392 if (ma!=nb || na!=mb) { 2393 *f = PETSC_FALSE; 2394 PetscFunctionReturn(0); 2395 } 2396 PetscCall(MatSeqAIJGetArrayRead(A,&va)); 2397 PetscCall(MatSeqAIJGetArrayRead(B,&vb)); 2398 aii = aij->i; bii = bij->i; 2399 adx = aij->j; bdx = bij->j; 2400 PetscCall(PetscMalloc1(ma,&aptr)); 2401 PetscCall(PetscMalloc1(mb,&bptr)); 2402 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2403 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2404 2405 *f = PETSC_TRUE; 2406 for (i=0; i<ma; i++) { 2407 while (aptr[i]<aii[i+1]) { 2408 PetscInt idc,idr; 2409 PetscScalar vc,vr; 2410 /* column/row index/value */ 2411 idc = adx[aptr[i]]; 2412 idr = bdx[bptr[idc]]; 2413 vc = va[aptr[i]]; 2414 vr = vb[bptr[idc]]; 2415 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2416 *f = PETSC_FALSE; 2417 goto done; 2418 } else { 2419 aptr[i]++; 2420 if (B || i!=idc) bptr[idc]++; 2421 } 2422 } 2423 } 2424 done: 2425 PetscCall(PetscFree(aptr)); 2426 PetscCall(PetscFree(bptr)); 2427 PetscCall(MatSeqAIJRestoreArrayRead(A,&va)); 2428 PetscCall(MatSeqAIJRestoreArrayRead(B,&vb)); 2429 PetscFunctionReturn(0); 2430 } 2431 2432 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2433 { 2434 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2435 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2436 MatScalar *va,*vb; 2437 PetscInt ma,na,mb,nb, i; 2438 2439 PetscFunctionBegin; 2440 PetscCall(MatGetSize(A,&ma,&na)); 2441 PetscCall(MatGetSize(B,&mb,&nb)); 2442 if (ma!=nb || na!=mb) { 2443 *f = PETSC_FALSE; 2444 PetscFunctionReturn(0); 2445 } 2446 aii = aij->i; bii = bij->i; 2447 adx = aij->j; bdx = bij->j; 2448 va = aij->a; vb = bij->a; 2449 PetscCall(PetscMalloc1(ma,&aptr)); 2450 PetscCall(PetscMalloc1(mb,&bptr)); 2451 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2452 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2453 2454 *f = PETSC_TRUE; 2455 for (i=0; i<ma; i++) { 2456 while (aptr[i]<aii[i+1]) { 2457 PetscInt idc,idr; 2458 PetscScalar vc,vr; 2459 /* column/row index/value */ 2460 idc = adx[aptr[i]]; 2461 idr = bdx[bptr[idc]]; 2462 vc = va[aptr[i]]; 2463 vr = vb[bptr[idc]]; 2464 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2465 *f = PETSC_FALSE; 2466 goto done; 2467 } else { 2468 aptr[i]++; 2469 if (B || i!=idc) bptr[idc]++; 2470 } 2471 } 2472 } 2473 done: 2474 PetscCall(PetscFree(aptr)); 2475 PetscCall(PetscFree(bptr)); 2476 PetscFunctionReturn(0); 2477 } 2478 2479 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2480 { 2481 PetscFunctionBegin; 2482 PetscCall(MatIsTranspose_SeqAIJ(A,A,tol,f)); 2483 PetscFunctionReturn(0); 2484 } 2485 2486 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2487 { 2488 PetscFunctionBegin; 2489 PetscCall(MatIsHermitianTranspose_SeqAIJ(A,A,tol,f)); 2490 PetscFunctionReturn(0); 2491 } 2492 2493 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2494 { 2495 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2496 const PetscScalar *l,*r; 2497 PetscScalar x; 2498 MatScalar *v; 2499 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2500 const PetscInt *jj; 2501 2502 PetscFunctionBegin; 2503 if (ll) { 2504 /* The local size is used so that VecMPI can be passed to this routine 2505 by MatDiagonalScale_MPIAIJ */ 2506 PetscCall(VecGetLocalSize(ll,&m)); 2507 PetscCheckFalse(m != A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2508 PetscCall(VecGetArrayRead(ll,&l)); 2509 PetscCall(MatSeqAIJGetArray(A,&v)); 2510 for (i=0; i<m; i++) { 2511 x = l[i]; 2512 M = a->i[i+1] - a->i[i]; 2513 for (j=0; j<M; j++) (*v++) *= x; 2514 } 2515 PetscCall(VecRestoreArrayRead(ll,&l)); 2516 PetscCall(PetscLogFlops(nz)); 2517 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2518 } 2519 if (rr) { 2520 PetscCall(VecGetLocalSize(rr,&n)); 2521 PetscCheckFalse(n != A->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2522 PetscCall(VecGetArrayRead(rr,&r)); 2523 PetscCall(MatSeqAIJGetArray(A,&v)); 2524 jj = a->j; 2525 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2526 PetscCall(MatSeqAIJRestoreArray(A,&v)); 2527 PetscCall(VecRestoreArrayRead(rr,&r)); 2528 PetscCall(PetscLogFlops(nz)); 2529 } 2530 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 2531 PetscFunctionReturn(0); 2532 } 2533 2534 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2535 { 2536 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2537 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2538 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2539 const PetscInt *irow,*icol; 2540 const PetscScalar *aa; 2541 PetscInt nrows,ncols; 2542 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2543 MatScalar *a_new,*mat_a; 2544 Mat C; 2545 PetscBool stride; 2546 2547 PetscFunctionBegin; 2548 PetscCall(ISGetIndices(isrow,&irow)); 2549 PetscCall(ISGetLocalSize(isrow,&nrows)); 2550 PetscCall(ISGetLocalSize(iscol,&ncols)); 2551 2552 PetscCall(PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride)); 2553 if (stride) { 2554 PetscCall(ISStrideGetInfo(iscol,&first,&step)); 2555 } else { 2556 first = 0; 2557 step = 0; 2558 } 2559 if (stride && step == 1) { 2560 /* special case of contiguous rows */ 2561 PetscCall(PetscMalloc2(nrows,&lens,nrows,&starts)); 2562 /* loop over new rows determining lens and starting points */ 2563 for (i=0; i<nrows; i++) { 2564 kstart = ai[irow[i]]; 2565 kend = kstart + ailen[irow[i]]; 2566 starts[i] = kstart; 2567 for (k=kstart; k<kend; k++) { 2568 if (aj[k] >= first) { 2569 starts[i] = k; 2570 break; 2571 } 2572 } 2573 sum = 0; 2574 while (k < kend) { 2575 if (aj[k++] >= first+ncols) break; 2576 sum++; 2577 } 2578 lens[i] = sum; 2579 } 2580 /* create submatrix */ 2581 if (scall == MAT_REUSE_MATRIX) { 2582 PetscInt n_cols,n_rows; 2583 PetscCall(MatGetSize(*B,&n_rows,&n_cols)); 2584 PetscCheckFalse(n_rows != nrows || n_cols != ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2585 PetscCall(MatZeroEntries(*B)); 2586 C = *B; 2587 } else { 2588 PetscInt rbs,cbs; 2589 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2590 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2591 PetscCall(ISGetBlockSize(isrow,&rbs)); 2592 PetscCall(ISGetBlockSize(iscol,&cbs)); 2593 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2594 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2595 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2596 } 2597 c = (Mat_SeqAIJ*)C->data; 2598 2599 /* loop over rows inserting into submatrix */ 2600 a_new = c->a; 2601 j_new = c->j; 2602 i_new = c->i; 2603 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2604 for (i=0; i<nrows; i++) { 2605 ii = starts[i]; 2606 lensi = lens[i]; 2607 for (k=0; k<lensi; k++) { 2608 *j_new++ = aj[ii+k] - first; 2609 } 2610 PetscCall(PetscArraycpy(a_new,aa + starts[i],lensi)); 2611 a_new += lensi; 2612 i_new[i+1] = i_new[i] + lensi; 2613 c->ilen[i] = lensi; 2614 } 2615 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2616 PetscCall(PetscFree2(lens,starts)); 2617 } else { 2618 PetscCall(ISGetIndices(iscol,&icol)); 2619 PetscCall(PetscCalloc1(oldcols,&smap)); 2620 PetscCall(PetscMalloc1(1+nrows,&lens)); 2621 for (i=0; i<ncols; i++) { 2622 PetscCheck(icol[i] < oldcols,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Requesting column beyond largest column icol[%" PetscInt_FMT "] %" PetscInt_FMT " >= A->cmap->n %" PetscInt_FMT,i,icol[i],oldcols); 2623 smap[icol[i]] = i+1; 2624 } 2625 2626 /* determine lens of each row */ 2627 for (i=0; i<nrows; i++) { 2628 kstart = ai[irow[i]]; 2629 kend = kstart + a->ilen[irow[i]]; 2630 lens[i] = 0; 2631 for (k=kstart; k<kend; k++) { 2632 if (smap[aj[k]]) { 2633 lens[i]++; 2634 } 2635 } 2636 } 2637 /* Create and fill new matrix */ 2638 if (scall == MAT_REUSE_MATRIX) { 2639 PetscBool equal; 2640 2641 c = (Mat_SeqAIJ*)((*B)->data); 2642 PetscCheckFalse((*B)->rmap->n != nrows || (*B)->cmap->n != ncols,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2643 PetscCall(PetscArraycmp(c->ilen,lens,(*B)->rmap->n,&equal)); 2644 PetscCheck(equal,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2645 PetscCall(PetscArrayzero(c->ilen,(*B)->rmap->n)); 2646 C = *B; 2647 } else { 2648 PetscInt rbs,cbs; 2649 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),&C)); 2650 PetscCall(MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE)); 2651 PetscCall(ISGetBlockSize(isrow,&rbs)); 2652 PetscCall(ISGetBlockSize(iscol,&cbs)); 2653 PetscCall(MatSetBlockSizes(C,rbs,cbs)); 2654 PetscCall(MatSetType(C,((PetscObject)A)->type_name)); 2655 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens)); 2656 } 2657 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 2658 c = (Mat_SeqAIJ*)(C->data); 2659 for (i=0; i<nrows; i++) { 2660 row = irow[i]; 2661 kstart = ai[row]; 2662 kend = kstart + a->ilen[row]; 2663 mat_i = c->i[i]; 2664 mat_j = c->j + mat_i; 2665 mat_a = c->a + mat_i; 2666 mat_ilen = c->ilen + i; 2667 for (k=kstart; k<kend; k++) { 2668 if ((tcol=smap[a->j[k]])) { 2669 *mat_j++ = tcol - 1; 2670 *mat_a++ = aa[k]; 2671 (*mat_ilen)++; 2672 2673 } 2674 } 2675 } 2676 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 2677 /* Free work space */ 2678 PetscCall(ISRestoreIndices(iscol,&icol)); 2679 PetscCall(PetscFree(smap)); 2680 PetscCall(PetscFree(lens)); 2681 /* sort */ 2682 for (i = 0; i < nrows; i++) { 2683 PetscInt ilen; 2684 2685 mat_i = c->i[i]; 2686 mat_j = c->j + mat_i; 2687 mat_a = c->a + mat_i; 2688 ilen = c->ilen[i]; 2689 PetscCall(PetscSortIntWithScalarArray(ilen,mat_j,mat_a)); 2690 } 2691 } 2692 #if defined(PETSC_HAVE_DEVICE) 2693 PetscCall(MatBindToCPU(C,A->boundtocpu)); 2694 #endif 2695 PetscCall(MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY)); 2696 PetscCall(MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY)); 2697 2698 PetscCall(ISRestoreIndices(isrow,&irow)); 2699 *B = C; 2700 PetscFunctionReturn(0); 2701 } 2702 2703 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2704 { 2705 Mat B; 2706 2707 PetscFunctionBegin; 2708 if (scall == MAT_INITIAL_MATRIX) { 2709 PetscCall(MatCreate(subComm,&B)); 2710 PetscCall(MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n)); 2711 PetscCall(MatSetBlockSizesFromMats(B,mat,mat)); 2712 PetscCall(MatSetType(B,MATSEQAIJ)); 2713 PetscCall(MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE)); 2714 *subMat = B; 2715 } else { 2716 PetscCall(MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN)); 2717 } 2718 PetscFunctionReturn(0); 2719 } 2720 2721 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2722 { 2723 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2724 Mat outA; 2725 PetscBool row_identity,col_identity; 2726 2727 PetscFunctionBegin; 2728 PetscCheckFalse(info->levels != 0,PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2729 2730 PetscCall(ISIdentity(row,&row_identity)); 2731 PetscCall(ISIdentity(col,&col_identity)); 2732 2733 outA = inA; 2734 outA->factortype = MAT_FACTOR_LU; 2735 PetscCall(PetscFree(inA->solvertype)); 2736 PetscCall(PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype)); 2737 2738 PetscCall(PetscObjectReference((PetscObject)row)); 2739 PetscCall(ISDestroy(&a->row)); 2740 2741 a->row = row; 2742 2743 PetscCall(PetscObjectReference((PetscObject)col)); 2744 PetscCall(ISDestroy(&a->col)); 2745 2746 a->col = col; 2747 2748 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2749 PetscCall(ISDestroy(&a->icol)); 2750 PetscCall(ISInvertPermutation(col,PETSC_DECIDE,&a->icol)); 2751 PetscCall(PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol)); 2752 2753 if (!a->solve_work) { /* this matrix may have been factored before */ 2754 PetscCall(PetscMalloc1(inA->rmap->n+1,&a->solve_work)); 2755 PetscCall(PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar))); 2756 } 2757 2758 PetscCall(MatMarkDiagonal_SeqAIJ(inA)); 2759 if (row_identity && col_identity) { 2760 PetscCall(MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info)); 2761 } else { 2762 PetscCall(MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info)); 2763 } 2764 PetscFunctionReturn(0); 2765 } 2766 2767 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2768 { 2769 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2770 PetscScalar *v; 2771 PetscBLASInt one = 1,bnz; 2772 2773 PetscFunctionBegin; 2774 PetscCall(MatSeqAIJGetArray(inA,&v)); 2775 PetscCall(PetscBLASIntCast(a->nz,&bnz)); 2776 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&alpha,v,&one)); 2777 PetscCall(PetscLogFlops(a->nz)); 2778 PetscCall(MatSeqAIJRestoreArray(inA,&v)); 2779 PetscCall(MatSeqAIJInvalidateDiagonal(inA)); 2780 PetscFunctionReturn(0); 2781 } 2782 2783 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2784 { 2785 PetscInt i; 2786 2787 PetscFunctionBegin; 2788 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2789 PetscCall(PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr)); 2790 2791 for (i=0; i<submatj->nrqr; ++i) { 2792 PetscCall(PetscFree(submatj->sbuf2[i])); 2793 } 2794 PetscCall(PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1)); 2795 2796 if (submatj->rbuf1) { 2797 PetscCall(PetscFree(submatj->rbuf1[0])); 2798 PetscCall(PetscFree(submatj->rbuf1)); 2799 } 2800 2801 for (i=0; i<submatj->nrqs; ++i) { 2802 PetscCall(PetscFree(submatj->rbuf3[i])); 2803 } 2804 PetscCall(PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3)); 2805 PetscCall(PetscFree(submatj->pa)); 2806 } 2807 2808 #if defined(PETSC_USE_CTABLE) 2809 PetscCall(PetscTableDestroy((PetscTable*)&submatj->rmap)); 2810 if (submatj->cmap_loc) PetscCall(PetscFree(submatj->cmap_loc)); 2811 PetscCall(PetscFree(submatj->rmap_loc)); 2812 #else 2813 PetscCall(PetscFree(submatj->rmap)); 2814 #endif 2815 2816 if (!submatj->allcolumns) { 2817 #if defined(PETSC_USE_CTABLE) 2818 PetscCall(PetscTableDestroy((PetscTable*)&submatj->cmap)); 2819 #else 2820 PetscCall(PetscFree(submatj->cmap)); 2821 #endif 2822 } 2823 PetscCall(PetscFree(submatj->row2proc)); 2824 2825 PetscCall(PetscFree(submatj)); 2826 PetscFunctionReturn(0); 2827 } 2828 2829 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2830 { 2831 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2832 Mat_SubSppt *submatj = c->submatis1; 2833 2834 PetscFunctionBegin; 2835 PetscCall((*submatj->destroy)(C)); 2836 PetscCall(MatDestroySubMatrix_Private(submatj)); 2837 PetscFunctionReturn(0); 2838 } 2839 2840 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2841 { 2842 PetscInt i; 2843 Mat C; 2844 Mat_SeqAIJ *c; 2845 Mat_SubSppt *submatj; 2846 2847 PetscFunctionBegin; 2848 for (i=0; i<n; i++) { 2849 C = (*mat)[i]; 2850 c = (Mat_SeqAIJ*)C->data; 2851 submatj = c->submatis1; 2852 if (submatj) { 2853 if (--((PetscObject)C)->refct <= 0) { 2854 PetscCall((*submatj->destroy)(C)); 2855 PetscCall(MatDestroySubMatrix_Private(submatj)); 2856 PetscCall(PetscFree(C->defaultvectype)); 2857 PetscCall(PetscLayoutDestroy(&C->rmap)); 2858 PetscCall(PetscLayoutDestroy(&C->cmap)); 2859 PetscCall(PetscHeaderDestroy(&C)); 2860 } 2861 } else { 2862 PetscCall(MatDestroy(&C)); 2863 } 2864 } 2865 2866 /* Destroy Dummy submatrices created for reuse */ 2867 PetscCall(MatDestroySubMatrices_Dummy(n,mat)); 2868 2869 PetscCall(PetscFree(*mat)); 2870 PetscFunctionReturn(0); 2871 } 2872 2873 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2874 { 2875 PetscInt i; 2876 2877 PetscFunctionBegin; 2878 if (scall == MAT_INITIAL_MATRIX) { 2879 PetscCall(PetscCalloc1(n+1,B)); 2880 } 2881 2882 for (i=0; i<n; i++) { 2883 PetscCall(MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i])); 2884 } 2885 PetscFunctionReturn(0); 2886 } 2887 2888 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2889 { 2890 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2891 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2892 const PetscInt *idx; 2893 PetscInt start,end,*ai,*aj; 2894 PetscBT table; 2895 2896 PetscFunctionBegin; 2897 m = A->rmap->n; 2898 ai = a->i; 2899 aj = a->j; 2900 2901 PetscCheckFalse(ov < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2902 2903 PetscCall(PetscMalloc1(m+1,&nidx)); 2904 PetscCall(PetscBTCreate(m,&table)); 2905 2906 for (i=0; i<is_max; i++) { 2907 /* Initialize the two local arrays */ 2908 isz = 0; 2909 PetscCall(PetscBTMemzero(m,table)); 2910 2911 /* Extract the indices, assume there can be duplicate entries */ 2912 PetscCall(ISGetIndices(is[i],&idx)); 2913 PetscCall(ISGetLocalSize(is[i],&n)); 2914 2915 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2916 for (j=0; j<n; ++j) { 2917 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2918 } 2919 PetscCall(ISRestoreIndices(is[i],&idx)); 2920 PetscCall(ISDestroy(&is[i])); 2921 2922 k = 0; 2923 for (j=0; j<ov; j++) { /* for each overlap */ 2924 n = isz; 2925 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2926 row = nidx[k]; 2927 start = ai[row]; 2928 end = ai[row+1]; 2929 for (l = start; l<end; l++) { 2930 val = aj[l]; 2931 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2932 } 2933 } 2934 } 2935 PetscCall(ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i))); 2936 } 2937 PetscCall(PetscBTDestroy(&table)); 2938 PetscCall(PetscFree(nidx)); 2939 PetscFunctionReturn(0); 2940 } 2941 2942 /* -------------------------------------------------------------- */ 2943 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2944 { 2945 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2946 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2947 const PetscInt *row,*col; 2948 PetscInt *cnew,j,*lens; 2949 IS icolp,irowp; 2950 PetscInt *cwork = NULL; 2951 PetscScalar *vwork = NULL; 2952 2953 PetscFunctionBegin; 2954 PetscCall(ISInvertPermutation(rowp,PETSC_DECIDE,&irowp)); 2955 PetscCall(ISGetIndices(irowp,&row)); 2956 PetscCall(ISInvertPermutation(colp,PETSC_DECIDE,&icolp)); 2957 PetscCall(ISGetIndices(icolp,&col)); 2958 2959 /* determine lengths of permuted rows */ 2960 PetscCall(PetscMalloc1(m+1,&lens)); 2961 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2962 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 2963 PetscCall(MatSetSizes(*B,m,n,m,n)); 2964 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 2965 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 2966 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens)); 2967 PetscCall(PetscFree(lens)); 2968 2969 PetscCall(PetscMalloc1(n,&cnew)); 2970 for (i=0; i<m; i++) { 2971 PetscCall(MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2972 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2973 PetscCall(MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES)); 2974 PetscCall(MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork)); 2975 } 2976 PetscCall(PetscFree(cnew)); 2977 2978 (*B)->assembled = PETSC_FALSE; 2979 2980 #if defined(PETSC_HAVE_DEVICE) 2981 PetscCall(MatBindToCPU(*B,A->boundtocpu)); 2982 #endif 2983 PetscCall(MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY)); 2984 PetscCall(MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY)); 2985 PetscCall(ISRestoreIndices(irowp,&row)); 2986 PetscCall(ISRestoreIndices(icolp,&col)); 2987 PetscCall(ISDestroy(&irowp)); 2988 PetscCall(ISDestroy(&icolp)); 2989 if (rowp == colp) { 2990 PetscCall(MatPropagateSymmetryOptions(A,*B)); 2991 } 2992 PetscFunctionReturn(0); 2993 } 2994 2995 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2996 { 2997 PetscFunctionBegin; 2998 /* If the two matrices have the same copy implementation, use fast copy. */ 2999 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 3000 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3001 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 3002 const PetscScalar *aa; 3003 3004 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 3005 PetscCheckFalse(a->i[A->rmap->n] != b->i[B->rmap->n],PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Number of nonzeros in two matrices are different %" PetscInt_FMT " != %" PetscInt_FMT,a->i[A->rmap->n],b->i[B->rmap->n]); 3006 PetscCall(PetscArraycpy(b->a,aa,a->i[A->rmap->n])); 3007 PetscCall(PetscObjectStateIncrease((PetscObject)B)); 3008 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 3009 } else { 3010 PetscCall(MatCopy_Basic(A,B,str)); 3011 } 3012 PetscFunctionReturn(0); 3013 } 3014 3015 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 3016 { 3017 PetscFunctionBegin; 3018 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,NULL)); 3019 PetscFunctionReturn(0); 3020 } 3021 3022 PETSC_INTERN PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 3023 { 3024 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3025 3026 PetscFunctionBegin; 3027 *array = a->a; 3028 PetscFunctionReturn(0); 3029 } 3030 3031 PETSC_INTERN PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 3032 { 3033 PetscFunctionBegin; 3034 *array = NULL; 3035 PetscFunctionReturn(0); 3036 } 3037 3038 /* 3039 Computes the number of nonzeros per row needed for preallocation when X and Y 3040 have different nonzero structure. 3041 */ 3042 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 3043 { 3044 PetscInt i,j,k,nzx,nzy; 3045 3046 PetscFunctionBegin; 3047 /* Set the number of nonzeros in the new matrix */ 3048 for (i=0; i<m; i++) { 3049 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 3050 nzx = xi[i+1] - xi[i]; 3051 nzy = yi[i+1] - yi[i]; 3052 nnz[i] = 0; 3053 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 3054 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 3055 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 3056 nnz[i]++; 3057 } 3058 for (; k<nzy; k++) nnz[i]++; 3059 } 3060 PetscFunctionReturn(0); 3061 } 3062 3063 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 3064 { 3065 PetscInt m = Y->rmap->N; 3066 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 3067 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 3068 3069 PetscFunctionBegin; 3070 /* Set the number of nonzeros in the new matrix */ 3071 PetscCall(MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz)); 3072 PetscFunctionReturn(0); 3073 } 3074 3075 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 3076 { 3077 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 3078 3079 PetscFunctionBegin; 3080 if (str == UNKNOWN_NONZERO_PATTERN || (PetscDefined(USE_DEBUG) && str == SAME_NONZERO_PATTERN)) { 3081 PetscBool e = x->nz == y->nz ? PETSC_TRUE : PETSC_FALSE; 3082 if (e) { 3083 PetscCall(PetscArraycmp(x->i,y->i,Y->rmap->n+1,&e)); 3084 if (e) { 3085 PetscCall(PetscArraycmp(x->j,y->j,y->nz,&e)); 3086 if (e) str = SAME_NONZERO_PATTERN; 3087 } 3088 } 3089 if (!e) PetscCheck(str != SAME_NONZERO_PATTERN,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"MatStructure is not SAME_NONZERO_PATTERN"); 3090 } 3091 if (str == SAME_NONZERO_PATTERN) { 3092 const PetscScalar *xa; 3093 PetscScalar *ya,alpha = a; 3094 PetscBLASInt one = 1,bnz; 3095 3096 PetscCall(PetscBLASIntCast(x->nz,&bnz)); 3097 PetscCall(MatSeqAIJGetArray(Y,&ya)); 3098 PetscCall(MatSeqAIJGetArrayRead(X,&xa)); 3099 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,xa,&one,ya,&one)); 3100 PetscCall(MatSeqAIJRestoreArrayRead(X,&xa)); 3101 PetscCall(MatSeqAIJRestoreArray(Y,&ya)); 3102 PetscCall(PetscLogFlops(2.0*bnz)); 3103 PetscCall(MatSeqAIJInvalidateDiagonal(Y)); 3104 PetscCall(PetscObjectStateIncrease((PetscObject)Y)); 3105 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 3106 PetscCall(MatAXPY_Basic(Y,a,X,str)); 3107 } else { 3108 Mat B; 3109 PetscInt *nnz; 3110 PetscCall(PetscMalloc1(Y->rmap->N,&nnz)); 3111 PetscCall(MatCreate(PetscObjectComm((PetscObject)Y),&B)); 3112 PetscCall(PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name)); 3113 PetscCall(MatSetLayouts(B,Y->rmap,Y->cmap)); 3114 PetscCall(MatSetType(B,((PetscObject)Y)->type_name)); 3115 PetscCall(MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz)); 3116 PetscCall(MatSeqAIJSetPreallocation(B,0,nnz)); 3117 PetscCall(MatAXPY_BasicWithPreallocation(B,Y,a,X,str)); 3118 PetscCall(MatHeaderMerge(Y,&B)); 3119 PetscCall(PetscFree(nnz)); 3120 } 3121 PetscFunctionReturn(0); 3122 } 3123 3124 PETSC_INTERN PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 3125 { 3126 #if defined(PETSC_USE_COMPLEX) 3127 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3128 PetscInt i,nz; 3129 PetscScalar *a; 3130 3131 PetscFunctionBegin; 3132 nz = aij->nz; 3133 PetscCall(MatSeqAIJGetArray(mat,&a)); 3134 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 3135 PetscCall(MatSeqAIJRestoreArray(mat,&a)); 3136 #else 3137 PetscFunctionBegin; 3138 #endif 3139 PetscFunctionReturn(0); 3140 } 3141 3142 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3143 { 3144 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3145 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3146 PetscReal atmp; 3147 PetscScalar *x; 3148 const MatScalar *aa,*av; 3149 3150 PetscFunctionBegin; 3151 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3152 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3153 aa = av; 3154 ai = a->i; 3155 aj = a->j; 3156 3157 PetscCall(VecSet(v,0.0)); 3158 PetscCall(VecGetArrayWrite(v,&x)); 3159 PetscCall(VecGetLocalSize(v,&n)); 3160 PetscCheckFalse(n != A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3161 for (i=0; i<m; i++) { 3162 ncols = ai[1] - ai[0]; ai++; 3163 for (j=0; j<ncols; j++) { 3164 atmp = PetscAbsScalar(*aa); 3165 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 3166 aa++; aj++; 3167 } 3168 } 3169 PetscCall(VecRestoreArrayWrite(v,&x)); 3170 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3171 PetscFunctionReturn(0); 3172 } 3173 3174 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3175 { 3176 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3177 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3178 PetscScalar *x; 3179 const MatScalar *aa,*av; 3180 3181 PetscFunctionBegin; 3182 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3183 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3184 aa = av; 3185 ai = a->i; 3186 aj = a->j; 3187 3188 PetscCall(VecSet(v,0.0)); 3189 PetscCall(VecGetArrayWrite(v,&x)); 3190 PetscCall(VecGetLocalSize(v,&n)); 3191 PetscCheckFalse(n != A->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3192 for (i=0; i<m; i++) { 3193 ncols = ai[1] - ai[0]; ai++; 3194 if (ncols == A->cmap->n) { /* row is dense */ 3195 x[i] = *aa; if (idx) idx[i] = 0; 3196 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 3197 x[i] = 0.0; 3198 if (idx) { 3199 for (j=0; j<ncols; j++) { /* find first implicit 0.0 in the row */ 3200 if (aj[j] > j) { 3201 idx[i] = j; 3202 break; 3203 } 3204 } 3205 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3206 if (j==ncols && j < A->cmap->n) idx[i] = j; 3207 } 3208 } 3209 for (j=0; j<ncols; j++) { 3210 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3211 aa++; aj++; 3212 } 3213 } 3214 PetscCall(VecRestoreArrayWrite(v,&x)); 3215 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3216 PetscFunctionReturn(0); 3217 } 3218 3219 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3220 { 3221 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3222 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 3223 PetscScalar *x; 3224 const MatScalar *aa,*av; 3225 3226 PetscFunctionBegin; 3227 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3228 aa = av; 3229 ai = a->i; 3230 aj = a->j; 3231 3232 PetscCall(VecSet(v,0.0)); 3233 PetscCall(VecGetArrayWrite(v,&x)); 3234 PetscCall(VecGetLocalSize(v,&n)); 3235 PetscCheckFalse(n != m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector, %" PetscInt_FMT " vs. %" PetscInt_FMT " rows", m, n); 3236 for (i=0; i<m; i++) { 3237 ncols = ai[1] - ai[0]; ai++; 3238 if (ncols == A->cmap->n) { /* row is dense */ 3239 x[i] = *aa; if (idx) idx[i] = 0; 3240 } else { /* row is sparse so already KNOW minimum is 0.0 or higher */ 3241 x[i] = 0.0; 3242 if (idx) { /* find first implicit 0.0 in the row */ 3243 for (j=0; j<ncols; j++) { 3244 if (aj[j] > j) { 3245 idx[i] = j; 3246 break; 3247 } 3248 } 3249 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3250 if (j==ncols && j < A->cmap->n) idx[i] = j; 3251 } 3252 } 3253 for (j=0; j<ncols; j++) { 3254 if (PetscAbsScalar(x[i]) > PetscAbsScalar(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3255 aa++; aj++; 3256 } 3257 } 3258 PetscCall(VecRestoreArrayWrite(v,&x)); 3259 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3260 PetscFunctionReturn(0); 3261 } 3262 3263 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 3264 { 3265 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 3266 PetscInt i,j,m = A->rmap->n,ncols,n; 3267 const PetscInt *ai,*aj; 3268 PetscScalar *x; 3269 const MatScalar *aa,*av; 3270 3271 PetscFunctionBegin; 3272 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3273 PetscCall(MatSeqAIJGetArrayRead(A,&av)); 3274 aa = av; 3275 ai = a->i; 3276 aj = a->j; 3277 3278 PetscCall(VecSet(v,0.0)); 3279 PetscCall(VecGetArrayWrite(v,&x)); 3280 PetscCall(VecGetLocalSize(v,&n)); 3281 PetscCheckFalse(n != m,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 3282 for (i=0; i<m; i++) { 3283 ncols = ai[1] - ai[0]; ai++; 3284 if (ncols == A->cmap->n) { /* row is dense */ 3285 x[i] = *aa; if (idx) idx[i] = 0; 3286 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 3287 x[i] = 0.0; 3288 if (idx) { /* find first implicit 0.0 in the row */ 3289 for (j=0; j<ncols; j++) { 3290 if (aj[j] > j) { 3291 idx[i] = j; 3292 break; 3293 } 3294 } 3295 /* in case first implicit 0.0 in the row occurs at ncols-th column */ 3296 if (j==ncols && j < A->cmap->n) idx[i] = j; 3297 } 3298 } 3299 for (j=0; j<ncols; j++) { 3300 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 3301 aa++; aj++; 3302 } 3303 } 3304 PetscCall(VecRestoreArrayWrite(v,&x)); 3305 PetscCall(MatSeqAIJRestoreArrayRead(A,&av)); 3306 PetscFunctionReturn(0); 3307 } 3308 3309 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 3310 { 3311 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 3312 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 3313 MatScalar *diag,work[25],*v_work; 3314 const PetscReal shift = 0.0; 3315 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 3316 3317 PetscFunctionBegin; 3318 allowzeropivot = PetscNot(A->erroriffailure); 3319 if (a->ibdiagvalid) { 3320 if (values) *values = a->ibdiag; 3321 PetscFunctionReturn(0); 3322 } 3323 PetscCall(MatMarkDiagonal_SeqAIJ(A)); 3324 if (!a->ibdiag) { 3325 PetscCall(PetscMalloc1(bs2*mbs,&a->ibdiag)); 3326 PetscCall(PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar))); 3327 } 3328 diag = a->ibdiag; 3329 if (values) *values = a->ibdiag; 3330 /* factor and invert each block */ 3331 switch (bs) { 3332 case 1: 3333 for (i=0; i<mbs; i++) { 3334 PetscCall(MatGetValues(A,1,&i,1,&i,diag+i)); 3335 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3336 if (allowzeropivot) { 3337 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3338 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3339 A->factorerror_zeropivot_row = i; 3340 PetscCall(PetscInfo(A,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON)); 3341 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_MAT_LU_ZRPVT,"Zero pivot, row %" PetscInt_FMT " pivot %g tolerance %g",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON); 3342 } 3343 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3344 } 3345 break; 3346 case 2: 3347 for (i=0; i<mbs; i++) { 3348 ij[0] = 2*i; ij[1] = 2*i + 1; 3349 PetscCall(MatGetValues(A,2,ij,2,ij,diag)); 3350 PetscCall(PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected)); 3351 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3352 PetscCall(PetscKernel_A_gets_transpose_A_2(diag)); 3353 diag += 4; 3354 } 3355 break; 3356 case 3: 3357 for (i=0; i<mbs; i++) { 3358 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3359 PetscCall(MatGetValues(A,3,ij,3,ij,diag)); 3360 PetscCall(PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected)); 3361 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3362 PetscCall(PetscKernel_A_gets_transpose_A_3(diag)); 3363 diag += 9; 3364 } 3365 break; 3366 case 4: 3367 for (i=0; i<mbs; i++) { 3368 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3369 PetscCall(MatGetValues(A,4,ij,4,ij,diag)); 3370 PetscCall(PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected)); 3371 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3372 PetscCall(PetscKernel_A_gets_transpose_A_4(diag)); 3373 diag += 16; 3374 } 3375 break; 3376 case 5: 3377 for (i=0; i<mbs; i++) { 3378 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3379 PetscCall(MatGetValues(A,5,ij,5,ij,diag)); 3380 PetscCall(PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected)); 3381 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3382 PetscCall(PetscKernel_A_gets_transpose_A_5(diag)); 3383 diag += 25; 3384 } 3385 break; 3386 case 6: 3387 for (i=0; i<mbs; i++) { 3388 ij[0] = 6*i; ij[1] = 6*i + 1; ij[2] = 6*i + 2; ij[3] = 6*i + 3; ij[4] = 6*i + 4; ij[5] = 6*i + 5; 3389 PetscCall(MatGetValues(A,6,ij,6,ij,diag)); 3390 PetscCall(PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected)); 3391 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3392 PetscCall(PetscKernel_A_gets_transpose_A_6(diag)); 3393 diag += 36; 3394 } 3395 break; 3396 case 7: 3397 for (i=0; i<mbs; i++) { 3398 ij[0] = 7*i; ij[1] = 7*i + 1; ij[2] = 7*i + 2; ij[3] = 7*i + 3; ij[4] = 7*i + 4; ij[5] = 7*i + 5; ij[5] = 7*i + 6; 3399 PetscCall(MatGetValues(A,7,ij,7,ij,diag)); 3400 PetscCall(PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected)); 3401 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3402 PetscCall(PetscKernel_A_gets_transpose_A_7(diag)); 3403 diag += 49; 3404 } 3405 break; 3406 default: 3407 PetscCall(PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ)); 3408 for (i=0; i<mbs; i++) { 3409 for (j=0; j<bs; j++) { 3410 IJ[j] = bs*i + j; 3411 } 3412 PetscCall(MatGetValues(A,bs,IJ,bs,IJ,diag)); 3413 PetscCall(PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected)); 3414 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3415 PetscCall(PetscKernel_A_gets_transpose_A_N(diag,bs)); 3416 diag += bs2; 3417 } 3418 PetscCall(PetscFree3(v_work,v_pivots,IJ)); 3419 } 3420 a->ibdiagvalid = PETSC_TRUE; 3421 PetscFunctionReturn(0); 3422 } 3423 3424 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3425 { 3426 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3427 PetscScalar a,*aa; 3428 PetscInt m,n,i,j,col; 3429 3430 PetscFunctionBegin; 3431 if (!x->assembled) { 3432 PetscCall(MatGetSize(x,&m,&n)); 3433 for (i=0; i<m; i++) { 3434 for (j=0; j<aij->imax[i]; j++) { 3435 PetscCall(PetscRandomGetValue(rctx,&a)); 3436 col = (PetscInt)(n*PetscRealPart(a)); 3437 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3438 } 3439 } 3440 } else { 3441 PetscCall(MatSeqAIJGetArrayWrite(x,&aa)); 3442 for (i=0; i<aij->nz; i++) PetscCall(PetscRandomGetValue(rctx,aa+i)); 3443 PetscCall(MatSeqAIJRestoreArrayWrite(x,&aa)); 3444 } 3445 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3446 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3447 PetscFunctionReturn(0); 3448 } 3449 3450 /* Like MatSetRandom_SeqAIJ, but do not set values on columns in range of [low, high) */ 3451 PetscErrorCode MatSetRandomSkipColumnRange_SeqAIJ_Private(Mat x,PetscInt low,PetscInt high,PetscRandom rctx) 3452 { 3453 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3454 PetscScalar a; 3455 PetscInt m,n,i,j,col,nskip; 3456 3457 PetscFunctionBegin; 3458 nskip = high - low; 3459 PetscCall(MatGetSize(x,&m,&n)); 3460 n -= nskip; /* shrink number of columns where nonzeros can be set */ 3461 for (i=0; i<m; i++) { 3462 for (j=0; j<aij->imax[i]; j++) { 3463 PetscCall(PetscRandomGetValue(rctx,&a)); 3464 col = (PetscInt)(n*PetscRealPart(a)); 3465 if (col >= low) col += nskip; /* shift col rightward to skip the hole */ 3466 PetscCall(MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES)); 3467 } 3468 } 3469 PetscCall(MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY)); 3470 PetscCall(MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY)); 3471 PetscFunctionReturn(0); 3472 } 3473 3474 /* -------------------------------------------------------------------*/ 3475 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3476 MatGetRow_SeqAIJ, 3477 MatRestoreRow_SeqAIJ, 3478 MatMult_SeqAIJ, 3479 /* 4*/ MatMultAdd_SeqAIJ, 3480 MatMultTranspose_SeqAIJ, 3481 MatMultTransposeAdd_SeqAIJ, 3482 NULL, 3483 NULL, 3484 NULL, 3485 /* 10*/ NULL, 3486 MatLUFactor_SeqAIJ, 3487 NULL, 3488 MatSOR_SeqAIJ, 3489 MatTranspose_SeqAIJ, 3490 /*1 5*/ MatGetInfo_SeqAIJ, 3491 MatEqual_SeqAIJ, 3492 MatGetDiagonal_SeqAIJ, 3493 MatDiagonalScale_SeqAIJ, 3494 MatNorm_SeqAIJ, 3495 /* 20*/ NULL, 3496 MatAssemblyEnd_SeqAIJ, 3497 MatSetOption_SeqAIJ, 3498 MatZeroEntries_SeqAIJ, 3499 /* 24*/ MatZeroRows_SeqAIJ, 3500 NULL, 3501 NULL, 3502 NULL, 3503 NULL, 3504 /* 29*/ MatSetUp_SeqAIJ, 3505 NULL, 3506 NULL, 3507 NULL, 3508 NULL, 3509 /* 34*/ MatDuplicate_SeqAIJ, 3510 NULL, 3511 NULL, 3512 MatILUFactor_SeqAIJ, 3513 NULL, 3514 /* 39*/ MatAXPY_SeqAIJ, 3515 MatCreateSubMatrices_SeqAIJ, 3516 MatIncreaseOverlap_SeqAIJ, 3517 MatGetValues_SeqAIJ, 3518 MatCopy_SeqAIJ, 3519 /* 44*/ MatGetRowMax_SeqAIJ, 3520 MatScale_SeqAIJ, 3521 MatShift_SeqAIJ, 3522 MatDiagonalSet_SeqAIJ, 3523 MatZeroRowsColumns_SeqAIJ, 3524 /* 49*/ MatSetRandom_SeqAIJ, 3525 MatGetRowIJ_SeqAIJ, 3526 MatRestoreRowIJ_SeqAIJ, 3527 MatGetColumnIJ_SeqAIJ, 3528 MatRestoreColumnIJ_SeqAIJ, 3529 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3530 NULL, 3531 NULL, 3532 MatPermute_SeqAIJ, 3533 NULL, 3534 /* 59*/ NULL, 3535 MatDestroy_SeqAIJ, 3536 MatView_SeqAIJ, 3537 NULL, 3538 NULL, 3539 /* 64*/ NULL, 3540 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3541 NULL, 3542 NULL, 3543 NULL, 3544 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3545 MatGetRowMinAbs_SeqAIJ, 3546 NULL, 3547 NULL, 3548 NULL, 3549 /* 74*/ NULL, 3550 MatFDColoringApply_AIJ, 3551 NULL, 3552 NULL, 3553 NULL, 3554 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3555 NULL, 3556 NULL, 3557 NULL, 3558 MatLoad_SeqAIJ, 3559 /* 84*/ MatIsSymmetric_SeqAIJ, 3560 MatIsHermitian_SeqAIJ, 3561 NULL, 3562 NULL, 3563 NULL, 3564 /* 89*/ NULL, 3565 NULL, 3566 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3567 NULL, 3568 NULL, 3569 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3570 NULL, 3571 NULL, 3572 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3573 NULL, 3574 /* 99*/ MatProductSetFromOptions_SeqAIJ, 3575 NULL, 3576 NULL, 3577 MatConjugate_SeqAIJ, 3578 NULL, 3579 /*104*/ MatSetValuesRow_SeqAIJ, 3580 MatRealPart_SeqAIJ, 3581 MatImaginaryPart_SeqAIJ, 3582 NULL, 3583 NULL, 3584 /*109*/ MatMatSolve_SeqAIJ, 3585 NULL, 3586 MatGetRowMin_SeqAIJ, 3587 NULL, 3588 MatMissingDiagonal_SeqAIJ, 3589 /*114*/ NULL, 3590 NULL, 3591 NULL, 3592 NULL, 3593 NULL, 3594 /*119*/ NULL, 3595 NULL, 3596 NULL, 3597 NULL, 3598 MatGetMultiProcBlock_SeqAIJ, 3599 /*124*/ MatFindNonzeroRows_SeqAIJ, 3600 MatGetColumnReductions_SeqAIJ, 3601 MatInvertBlockDiagonal_SeqAIJ, 3602 MatInvertVariableBlockDiagonal_SeqAIJ, 3603 NULL, 3604 /*129*/ NULL, 3605 NULL, 3606 NULL, 3607 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3608 MatTransposeColoringCreate_SeqAIJ, 3609 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3610 MatTransColoringApplyDenToSp_SeqAIJ, 3611 NULL, 3612 NULL, 3613 MatRARtNumeric_SeqAIJ_SeqAIJ, 3614 /*139*/NULL, 3615 NULL, 3616 NULL, 3617 MatFDColoringSetUp_SeqXAIJ, 3618 MatFindOffBlockDiagonalEntries_SeqAIJ, 3619 MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3620 /*145*/MatDestroySubMatrices_SeqAIJ, 3621 NULL, 3622 NULL 3623 }; 3624 3625 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3626 { 3627 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3628 PetscInt i,nz,n; 3629 3630 PetscFunctionBegin; 3631 nz = aij->maxnz; 3632 n = mat->rmap->n; 3633 for (i=0; i<nz; i++) { 3634 aij->j[i] = indices[i]; 3635 } 3636 aij->nz = nz; 3637 for (i=0; i<n; i++) { 3638 aij->ilen[i] = aij->imax[i]; 3639 } 3640 PetscFunctionReturn(0); 3641 } 3642 3643 /* 3644 * Given a sparse matrix with global column indices, compact it by using a local column space. 3645 * The result matrix helps saving memory in other algorithms, such as MatPtAPSymbolic_MPIAIJ_MPIAIJ_scalable() 3646 */ 3647 PetscErrorCode MatSeqAIJCompactOutExtraColumns_SeqAIJ(Mat mat, ISLocalToGlobalMapping *mapping) 3648 { 3649 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3650 PetscTable gid1_lid1; 3651 PetscTablePosition tpos; 3652 PetscInt gid,lid,i,ec,nz = aij->nz; 3653 PetscInt *garray,*jj = aij->j; 3654 3655 PetscFunctionBegin; 3656 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3657 PetscValidPointer(mapping,2); 3658 /* use a table */ 3659 PetscCall(PetscTableCreate(mat->rmap->n,mat->cmap->N+1,&gid1_lid1)); 3660 ec = 0; 3661 for (i=0; i<nz; i++) { 3662 PetscInt data,gid1 = jj[i] + 1; 3663 PetscCall(PetscTableFind(gid1_lid1,gid1,&data)); 3664 if (!data) { 3665 /* one based table */ 3666 PetscCall(PetscTableAdd(gid1_lid1,gid1,++ec,INSERT_VALUES)); 3667 } 3668 } 3669 /* form array of columns we need */ 3670 PetscCall(PetscMalloc1(ec,&garray)); 3671 PetscCall(PetscTableGetHeadPosition(gid1_lid1,&tpos)); 3672 while (tpos) { 3673 PetscCall(PetscTableGetNext(gid1_lid1,&tpos,&gid,&lid)); 3674 gid--; 3675 lid--; 3676 garray[lid] = gid; 3677 } 3678 PetscCall(PetscSortInt(ec,garray)); /* sort, and rebuild */ 3679 PetscCall(PetscTableRemoveAll(gid1_lid1)); 3680 for (i=0; i<ec; i++) { 3681 PetscCall(PetscTableAdd(gid1_lid1,garray[i]+1,i+1,INSERT_VALUES)); 3682 } 3683 /* compact out the extra columns in B */ 3684 for (i=0; i<nz; i++) { 3685 PetscInt gid1 = jj[i] + 1; 3686 PetscCall(PetscTableFind(gid1_lid1,gid1,&lid)); 3687 lid--; 3688 jj[i] = lid; 3689 } 3690 PetscCall(PetscLayoutDestroy(&mat->cmap)); 3691 PetscCall(PetscTableDestroy(&gid1_lid1)); 3692 PetscCall(PetscLayoutCreateFromSizes(PetscObjectComm((PetscObject)mat),ec,ec,1,&mat->cmap)); 3693 PetscCall(ISLocalToGlobalMappingCreate(PETSC_COMM_SELF,mat->cmap->bs,mat->cmap->n,garray,PETSC_OWN_POINTER,mapping)); 3694 PetscCall(ISLocalToGlobalMappingSetType(*mapping,ISLOCALTOGLOBALMAPPINGHASH)); 3695 PetscFunctionReturn(0); 3696 } 3697 3698 /*@ 3699 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3700 in the matrix. 3701 3702 Input Parameters: 3703 + mat - the SeqAIJ matrix 3704 - indices - the column indices 3705 3706 Level: advanced 3707 3708 Notes: 3709 This can be called if you have precomputed the nonzero structure of the 3710 matrix and want to provide it to the matrix object to improve the performance 3711 of the MatSetValues() operation. 3712 3713 You MUST have set the correct numbers of nonzeros per row in the call to 3714 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3715 3716 MUST be called before any calls to MatSetValues(); 3717 3718 The indices should start with zero, not one. 3719 3720 @*/ 3721 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3722 { 3723 PetscFunctionBegin; 3724 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3725 PetscValidIntPointer(indices,2); 3726 PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices)); 3727 PetscFunctionReturn(0); 3728 } 3729 3730 /* ----------------------------------------------------------------------------------------*/ 3731 3732 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3733 { 3734 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3735 size_t nz = aij->i[mat->rmap->n]; 3736 3737 PetscFunctionBegin; 3738 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3739 3740 /* allocate space for values if not already there */ 3741 if (!aij->saved_values) { 3742 PetscCall(PetscMalloc1(nz+1,&aij->saved_values)); 3743 PetscCall(PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar))); 3744 } 3745 3746 /* copy values over */ 3747 PetscCall(PetscArraycpy(aij->saved_values,aij->a,nz)); 3748 PetscFunctionReturn(0); 3749 } 3750 3751 /*@ 3752 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3753 example, reuse of the linear part of a Jacobian, while recomputing the 3754 nonlinear portion. 3755 3756 Collect on Mat 3757 3758 Input Parameters: 3759 . mat - the matrix (currently only AIJ matrices support this option) 3760 3761 Level: advanced 3762 3763 Common Usage, with SNESSolve(): 3764 $ Create Jacobian matrix 3765 $ Set linear terms into matrix 3766 $ Apply boundary conditions to matrix, at this time matrix must have 3767 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3768 $ boundary conditions again will not change the nonzero structure 3769 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3770 $ ierr = MatStoreValues(mat); 3771 $ Call SNESSetJacobian() with matrix 3772 $ In your Jacobian routine 3773 $ ierr = MatRetrieveValues(mat); 3774 $ Set nonlinear terms in matrix 3775 3776 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3777 $ // build linear portion of Jacobian 3778 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3779 $ ierr = MatStoreValues(mat); 3780 $ loop over nonlinear iterations 3781 $ ierr = MatRetrieveValues(mat); 3782 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3783 $ // call MatAssemblyBegin/End() on matrix 3784 $ Solve linear system with Jacobian 3785 $ endloop 3786 3787 Notes: 3788 Matrix must already be assemblied before calling this routine 3789 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3790 calling this routine. 3791 3792 When this is called multiple times it overwrites the previous set of stored values 3793 and does not allocated additional space. 3794 3795 .seealso: MatRetrieveValues() 3796 3797 @*/ 3798 PetscErrorCode MatStoreValues(Mat mat) 3799 { 3800 PetscFunctionBegin; 3801 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3802 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3803 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3804 PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat)); 3805 PetscFunctionReturn(0); 3806 } 3807 3808 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3809 { 3810 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3811 PetscInt nz = aij->i[mat->rmap->n]; 3812 3813 PetscFunctionBegin; 3814 PetscCheck(aij->nonew,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3815 PetscCheck(aij->saved_values,PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3816 /* copy values over */ 3817 PetscCall(PetscArraycpy(aij->a,aij->saved_values,nz)); 3818 PetscFunctionReturn(0); 3819 } 3820 3821 /*@ 3822 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3823 example, reuse of the linear part of a Jacobian, while recomputing the 3824 nonlinear portion. 3825 3826 Collect on Mat 3827 3828 Input Parameters: 3829 . mat - the matrix (currently only AIJ matrices support this option) 3830 3831 Level: advanced 3832 3833 .seealso: MatStoreValues() 3834 3835 @*/ 3836 PetscErrorCode MatRetrieveValues(Mat mat) 3837 { 3838 PetscFunctionBegin; 3839 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3840 PetscCheck(mat->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3841 PetscCheck(!mat->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3842 PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat)); 3843 PetscFunctionReturn(0); 3844 } 3845 3846 /* --------------------------------------------------------------------------------*/ 3847 /*@C 3848 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3849 (the default parallel PETSc format). For good matrix assembly performance 3850 the user should preallocate the matrix storage by setting the parameter nz 3851 (or the array nnz). By setting these parameters accurately, performance 3852 during matrix assembly can be increased by more than a factor of 50. 3853 3854 Collective 3855 3856 Input Parameters: 3857 + comm - MPI communicator, set to PETSC_COMM_SELF 3858 . m - number of rows 3859 . n - number of columns 3860 . nz - number of nonzeros per row (same for all rows) 3861 - nnz - array containing the number of nonzeros in the various rows 3862 (possibly different for each row) or NULL 3863 3864 Output Parameter: 3865 . A - the matrix 3866 3867 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3868 MatXXXXSetPreallocation() paradigm instead of this routine directly. 3869 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3870 3871 Notes: 3872 If nnz is given then nz is ignored 3873 3874 The AIJ format (also called the Yale sparse matrix format or 3875 compressed row storage), is fully compatible with standard Fortran 77 3876 storage. That is, the stored row and column indices can begin at 3877 either one (as in Fortran) or zero. See the users' manual for details. 3878 3879 Specify the preallocated storage with either nz or nnz (not both). 3880 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3881 allocation. For large problems you MUST preallocate memory or you 3882 will get TERRIBLE performance, see the users' manual chapter on matrices. 3883 3884 By default, this format uses inodes (identical nodes) when possible, to 3885 improve numerical efficiency of matrix-vector products and solves. We 3886 search for consecutive rows with the same nonzero structure, thereby 3887 reusing matrix information to achieve increased efficiency. 3888 3889 Options Database Keys: 3890 + -mat_no_inode - Do not use inodes 3891 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3892 3893 Level: intermediate 3894 3895 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3896 3897 @*/ 3898 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3899 { 3900 PetscFunctionBegin; 3901 PetscCall(MatCreate(comm,A)); 3902 PetscCall(MatSetSizes(*A,m,n,m,n)); 3903 PetscCall(MatSetType(*A,MATSEQAIJ)); 3904 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz)); 3905 PetscFunctionReturn(0); 3906 } 3907 3908 /*@C 3909 MatSeqAIJSetPreallocation - For good matrix assembly performance 3910 the user should preallocate the matrix storage by setting the parameter nz 3911 (or the array nnz). By setting these parameters accurately, performance 3912 during matrix assembly can be increased by more than a factor of 50. 3913 3914 Collective 3915 3916 Input Parameters: 3917 + B - The matrix 3918 . nz - number of nonzeros per row (same for all rows) 3919 - nnz - array containing the number of nonzeros in the various rows 3920 (possibly different for each row) or NULL 3921 3922 Notes: 3923 If nnz is given then nz is ignored 3924 3925 The AIJ format (also called the Yale sparse matrix format or 3926 compressed row storage), is fully compatible with standard Fortran 77 3927 storage. That is, the stored row and column indices can begin at 3928 either one (as in Fortran) or zero. See the users' manual for details. 3929 3930 Specify the preallocated storage with either nz or nnz (not both). 3931 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3932 allocation. For large problems you MUST preallocate memory or you 3933 will get TERRIBLE performance, see the users' manual chapter on matrices. 3934 3935 You can call MatGetInfo() to get information on how effective the preallocation was; 3936 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3937 You can also run with the option -info and look for messages with the string 3938 malloc in them to see if additional memory allocation was needed. 3939 3940 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3941 entries or columns indices 3942 3943 By default, this format uses inodes (identical nodes) when possible, to 3944 improve numerical efficiency of matrix-vector products and solves. We 3945 search for consecutive rows with the same nonzero structure, thereby 3946 reusing matrix information to achieve increased efficiency. 3947 3948 Options Database Keys: 3949 + -mat_no_inode - Do not use inodes 3950 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3951 3952 Level: intermediate 3953 3954 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo(), 3955 MatSeqAIJSetTotalPreallocation() 3956 3957 @*/ 3958 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3959 { 3960 PetscFunctionBegin; 3961 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3962 PetscValidType(B,1); 3963 PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz)); 3964 PetscFunctionReturn(0); 3965 } 3966 3967 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3968 { 3969 Mat_SeqAIJ *b; 3970 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3971 PetscInt i; 3972 3973 PetscFunctionBegin; 3974 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3975 if (nz == MAT_SKIP_ALLOCATION) { 3976 skipallocation = PETSC_TRUE; 3977 nz = 0; 3978 } 3979 PetscCall(PetscLayoutSetUp(B->rmap)); 3980 PetscCall(PetscLayoutSetUp(B->cmap)); 3981 3982 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3983 PetscCheckFalse(nz < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %" PetscInt_FMT,nz); 3984 if (PetscUnlikelyDebug(nnz)) { 3985 for (i=0; i<B->rmap->n; i++) { 3986 PetscCheckFalse(nnz[i] < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be less than 0: local row %" PetscInt_FMT " value %" PetscInt_FMT,i,nnz[i]); 3987 PetscCheckFalse(nnz[i] > B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nnz cannot be greater than row length: local row %" PetscInt_FMT " value %" PetscInt_FMT " rowlength %" PetscInt_FMT,i,nnz[i],B->cmap->n); 3988 } 3989 } 3990 3991 B->preallocated = PETSC_TRUE; 3992 3993 b = (Mat_SeqAIJ*)B->data; 3994 3995 if (!skipallocation) { 3996 if (!b->imax) { 3997 PetscCall(PetscMalloc1(B->rmap->n,&b->imax)); 3998 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 3999 } 4000 if (!b->ilen) { 4001 /* b->ilen will count nonzeros in each row so far. */ 4002 PetscCall(PetscCalloc1(B->rmap->n,&b->ilen)); 4003 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4004 } else { 4005 PetscCall(PetscMemzero(b->ilen,B->rmap->n*sizeof(PetscInt))); 4006 } 4007 if (!b->ipre) { 4008 PetscCall(PetscMalloc1(B->rmap->n,&b->ipre)); 4009 PetscCall(PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt))); 4010 } 4011 if (!nnz) { 4012 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 4013 else if (nz < 0) nz = 1; 4014 nz = PetscMin(nz,B->cmap->n); 4015 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 4016 nz = nz*B->rmap->n; 4017 } else { 4018 PetscInt64 nz64 = 0; 4019 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz64 += nnz[i];} 4020 PetscCall(PetscIntCast(nz64,&nz)); 4021 } 4022 4023 /* allocate the matrix space */ 4024 /* FIXME: should B's old memory be unlogged? */ 4025 PetscCall(MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i)); 4026 if (B->structure_only) { 4027 PetscCall(PetscMalloc1(nz,&b->j)); 4028 PetscCall(PetscMalloc1(B->rmap->n+1,&b->i)); 4029 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt))); 4030 } else { 4031 PetscCall(PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i)); 4032 PetscCall(PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)))); 4033 } 4034 b->i[0] = 0; 4035 for (i=1; i<B->rmap->n+1; i++) { 4036 b->i[i] = b->i[i-1] + b->imax[i-1]; 4037 } 4038 if (B->structure_only) { 4039 b->singlemalloc = PETSC_FALSE; 4040 b->free_a = PETSC_FALSE; 4041 } else { 4042 b->singlemalloc = PETSC_TRUE; 4043 b->free_a = PETSC_TRUE; 4044 } 4045 b->free_ij = PETSC_TRUE; 4046 } else { 4047 b->free_a = PETSC_FALSE; 4048 b->free_ij = PETSC_FALSE; 4049 } 4050 4051 if (b->ipre && nnz != b->ipre && b->imax) { 4052 /* reserve user-requested sparsity */ 4053 PetscCall(PetscArraycpy(b->ipre,b->imax,B->rmap->n)); 4054 } 4055 4056 b->nz = 0; 4057 b->maxnz = nz; 4058 B->info.nz_unneeded = (double)b->maxnz; 4059 if (realalloc) { 4060 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE)); 4061 } 4062 B->was_assembled = PETSC_FALSE; 4063 B->assembled = PETSC_FALSE; 4064 PetscFunctionReturn(0); 4065 } 4066 4067 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 4068 { 4069 Mat_SeqAIJ *a; 4070 PetscInt i; 4071 4072 PetscFunctionBegin; 4073 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4074 4075 /* Check local size. If zero, then return */ 4076 if (!A->rmap->n) PetscFunctionReturn(0); 4077 4078 a = (Mat_SeqAIJ*)A->data; 4079 /* if no saved info, we error out */ 4080 PetscCheck(a->ipre,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"No saved preallocation info "); 4081 4082 PetscCheckFalse(!a->i || !a->j || !a->a || !a->imax || !a->ilen,PETSC_COMM_SELF,PETSC_ERR_ARG_NULL,"Memory info is incomplete, and can not reset preallocation "); 4083 4084 PetscCall(PetscArraycpy(a->imax,a->ipre,A->rmap->n)); 4085 PetscCall(PetscArrayzero(a->ilen,A->rmap->n)); 4086 a->i[0] = 0; 4087 for (i=1; i<A->rmap->n+1; i++) { 4088 a->i[i] = a->i[i-1] + a->imax[i-1]; 4089 } 4090 A->preallocated = PETSC_TRUE; 4091 a->nz = 0; 4092 a->maxnz = a->i[A->rmap->n]; 4093 A->info.nz_unneeded = (double)a->maxnz; 4094 A->was_assembled = PETSC_FALSE; 4095 A->assembled = PETSC_FALSE; 4096 PetscFunctionReturn(0); 4097 } 4098 4099 /*@ 4100 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 4101 4102 Input Parameters: 4103 + B - the matrix 4104 . i - the indices into j for the start of each row (starts with zero) 4105 . j - the column indices for each row (starts with zero) these must be sorted for each row 4106 - v - optional values in the matrix 4107 4108 Level: developer 4109 4110 Notes: 4111 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 4112 4113 This routine may be called multiple times with different nonzero patterns (or the same nonzero pattern). The nonzero 4114 structure will be the union of all the previous nonzero structures. 4115 4116 Developer Notes: 4117 An optimization could be added to the implementation where it checks if the i, and j are identical to the current i and j and 4118 then just copies the v values directly with PetscMemcpy(). 4119 4120 This routine could also take a PetscCopyMode argument to allow sharing the values instead of always copying them. 4121 4122 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ, MatResetPreallocation() 4123 @*/ 4124 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 4125 { 4126 PetscFunctionBegin; 4127 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 4128 PetscValidType(B,1); 4129 PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v)); 4130 PetscFunctionReturn(0); 4131 } 4132 4133 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 4134 { 4135 PetscInt i; 4136 PetscInt m,n; 4137 PetscInt nz; 4138 PetscInt *nnz; 4139 4140 PetscFunctionBegin; 4141 PetscCheckFalse(Ii[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %" PetscInt_FMT, Ii[0]); 4142 4143 PetscCall(PetscLayoutSetUp(B->rmap)); 4144 PetscCall(PetscLayoutSetUp(B->cmap)); 4145 4146 PetscCall(MatGetSize(B, &m, &n)); 4147 PetscCall(PetscMalloc1(m+1, &nnz)); 4148 for (i = 0; i < m; i++) { 4149 nz = Ii[i+1]- Ii[i]; 4150 PetscCheckFalse(nz < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %" PetscInt_FMT " has a negative number of columns %" PetscInt_FMT, i, nz); 4151 nnz[i] = nz; 4152 } 4153 PetscCall(MatSeqAIJSetPreallocation(B, 0, nnz)); 4154 PetscCall(PetscFree(nnz)); 4155 4156 for (i = 0; i < m; i++) { 4157 PetscCall(MatSetValues_SeqAIJ(B, 1, &i, Ii[i+1] - Ii[i], J+Ii[i], v ? v + Ii[i] : NULL, INSERT_VALUES)); 4158 } 4159 4160 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 4161 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 4162 4163 PetscCall(MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 4164 PetscFunctionReturn(0); 4165 } 4166 4167 /*@ 4168 MatSeqAIJKron - Computes C, the Kronecker product of A and B. 4169 4170 Input Parameters: 4171 + A - left-hand side matrix 4172 . B - right-hand side matrix 4173 - reuse - either MAT_INITIAL_MATRIX or MAT_REUSE_MATRIX 4174 4175 Output Parameter: 4176 . C - Kronecker product of A and B 4177 4178 Level: intermediate 4179 4180 Notes: 4181 MAT_REUSE_MATRIX can only be used when the nonzero structure of the product matrix has not changed from that last call to MatSeqAIJKron(). 4182 4183 .seealso: MatCreateSeqAIJ(), MATSEQAIJ, MATKAIJ, MatReuse 4184 @*/ 4185 PetscErrorCode MatSeqAIJKron(Mat A,Mat B,MatReuse reuse,Mat *C) 4186 { 4187 PetscFunctionBegin; 4188 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 4189 PetscValidType(A,1); 4190 PetscValidHeaderSpecific(B,MAT_CLASSID,2); 4191 PetscValidType(B,2); 4192 PetscValidPointer(C,4); 4193 if (reuse == MAT_REUSE_MATRIX) { 4194 PetscValidHeaderSpecific(*C,MAT_CLASSID,4); 4195 PetscValidType(*C,4); 4196 } 4197 PetscTryMethod(A,"MatSeqAIJKron_C",(Mat,Mat,MatReuse,Mat*),(A,B,reuse,C)); 4198 PetscFunctionReturn(0); 4199 } 4200 4201 PetscErrorCode MatSeqAIJKron_SeqAIJ(Mat A,Mat B,MatReuse reuse,Mat *C) 4202 { 4203 Mat newmat; 4204 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4205 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 4206 PetscScalar *v; 4207 const PetscScalar *aa,*ba; 4208 PetscInt *i,*j,m,n,p,q,nnz = 0,am = A->rmap->n,bm = B->rmap->n,an = A->cmap->n, bn = B->cmap->n; 4209 PetscBool flg; 4210 4211 PetscFunctionBegin; 4212 PetscCheck(!A->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4213 PetscCheck(A->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4214 PetscCheck(!B->factortype,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 4215 PetscCheck(B->assembled,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 4216 PetscCall(PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&flg)); 4217 PetscCheck(flg,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatType %s",((PetscObject)B)->type_name); 4218 PetscCheckFalse(reuse != MAT_INITIAL_MATRIX && reuse != MAT_REUSE_MATRIX,PETSC_COMM_SELF,PETSC_ERR_SUP,"MatReuse %d",(int)reuse); 4219 if (reuse == MAT_INITIAL_MATRIX) { 4220 PetscCall(PetscMalloc2(am*bm+1,&i,a->i[am]*b->i[bm],&j)); 4221 PetscCall(MatCreate(PETSC_COMM_SELF,&newmat)); 4222 PetscCall(MatSetSizes(newmat,am*bm,an*bn,am*bm,an*bn)); 4223 PetscCall(MatSetType(newmat,MATAIJ)); 4224 i[0] = 0; 4225 for (m = 0; m < am; ++m) { 4226 for (p = 0; p < bm; ++p) { 4227 i[m*bm + p + 1] = i[m*bm + p] + (a->i[m+1] - a->i[m]) * (b->i[p+1] - b->i[p]); 4228 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4229 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4230 j[nnz++] = a->j[n]*bn + b->j[q]; 4231 } 4232 } 4233 } 4234 } 4235 PetscCall(MatSeqAIJSetPreallocationCSR(newmat,i,j,NULL)); 4236 *C = newmat; 4237 PetscCall(PetscFree2(i,j)); 4238 nnz = 0; 4239 } 4240 PetscCall(MatSeqAIJGetArray(*C,&v)); 4241 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4242 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 4243 for (m = 0; m < am; ++m) { 4244 for (p = 0; p < bm; ++p) { 4245 for (n = a->i[m]; n < a->i[m+1]; ++n) { 4246 for (q = b->i[p]; q < b->i[p+1]; ++q) { 4247 v[nnz++] = aa[n] * ba[q]; 4248 } 4249 } 4250 } 4251 } 4252 PetscCall(MatSeqAIJRestoreArray(*C,&v)); 4253 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 4254 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 4255 PetscFunctionReturn(0); 4256 } 4257 4258 #include <../src/mat/impls/dense/seq/dense.h> 4259 #include <petsc/private/kernels/petscaxpy.h> 4260 4261 /* 4262 Computes (B'*A')' since computing B*A directly is untenable 4263 4264 n p p 4265 [ ] [ ] [ ] 4266 m [ A ] * n [ B ] = m [ C ] 4267 [ ] [ ] [ ] 4268 4269 */ 4270 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 4271 { 4272 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 4273 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 4274 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 4275 PetscInt i,j,n,m,q,p; 4276 const PetscInt *ii,*idx; 4277 const PetscScalar *b,*a,*a_q; 4278 PetscScalar *c,*c_q; 4279 PetscInt clda = sub_c->lda; 4280 PetscInt alda = sub_a->lda; 4281 4282 PetscFunctionBegin; 4283 m = A->rmap->n; 4284 n = A->cmap->n; 4285 p = B->cmap->n; 4286 a = sub_a->v; 4287 b = sub_b->a; 4288 c = sub_c->v; 4289 if (clda == m) { 4290 PetscCall(PetscArrayzero(c,m*p)); 4291 } else { 4292 for (j=0;j<p;j++) 4293 for (i=0;i<m;i++) 4294 c[j*clda + i] = 0.0; 4295 } 4296 ii = sub_b->i; 4297 idx = sub_b->j; 4298 for (i=0; i<n; i++) { 4299 q = ii[i+1] - ii[i]; 4300 while (q-->0) { 4301 c_q = c + clda*(*idx); 4302 a_q = a + alda*i; 4303 PetscKernelAXPY(c_q,*b,a_q,m); 4304 idx++; 4305 b++; 4306 } 4307 } 4308 PetscFunctionReturn(0); 4309 } 4310 4311 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat C) 4312 { 4313 PetscInt m=A->rmap->n,n=B->cmap->n; 4314 PetscBool cisdense; 4315 4316 PetscFunctionBegin; 4317 PetscCheckFalse(A->cmap->n != B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"A->cmap->n %" PetscInt_FMT " != B->rmap->n %" PetscInt_FMT,A->cmap->n,B->rmap->n); 4318 PetscCall(MatSetSizes(C,m,n,m,n)); 4319 PetscCall(MatSetBlockSizesFromMats(C,A,B)); 4320 PetscCall(PetscObjectTypeCompareAny((PetscObject)C,&cisdense,MATSEQDENSE,MATSEQDENSECUDA,"")); 4321 if (!cisdense) { 4322 PetscCall(MatSetType(C,MATDENSE)); 4323 } 4324 PetscCall(MatSetUp(C)); 4325 4326 C->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 4327 PetscFunctionReturn(0); 4328 } 4329 4330 /* ----------------------------------------------------------------*/ 4331 /*MC 4332 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 4333 based on compressed sparse row format. 4334 4335 Options Database Keys: 4336 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 4337 4338 Level: beginner 4339 4340 Notes: 4341 MatSetValues() may be called for this matrix type with a NULL argument for the numerical values, 4342 in this case the values associated with the rows and columns one passes in are set to zero 4343 in the matrix 4344 4345 MatSetOptions(,MAT_STRUCTURE_ONLY,PETSC_TRUE) may be called for this matrix type. In this no 4346 space is allocated for the nonzero entries and any entries passed with MatSetValues() are ignored 4347 4348 Developer Notes: 4349 It would be nice if all matrix formats supported passing NULL in for the numerical values 4350 4351 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType, MATSELL, MATSEQSELL, MATMPISELL 4352 M*/ 4353 4354 /*MC 4355 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 4356 4357 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 4358 and MATMPIAIJ otherwise. As a result, for single process communicators, 4359 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4360 for communicators controlling multiple processes. It is recommended that you call both of 4361 the above preallocation routines for simplicity. 4362 4363 Options Database Keys: 4364 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 4365 4366 Developer Notes: 4367 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJSELL, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 4368 enough exist. 4369 4370 Level: beginner 4371 4372 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ, MATMPIAIJ, MATSELL, MATSEQSELL, MATMPISELL 4373 M*/ 4374 4375 /*MC 4376 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 4377 4378 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 4379 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 4380 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 4381 for communicators controlling multiple processes. It is recommended that you call both of 4382 the above preallocation routines for simplicity. 4383 4384 Options Database Keys: 4385 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 4386 4387 Level: beginner 4388 4389 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 4390 M*/ 4391 4392 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 4393 #if defined(PETSC_HAVE_ELEMENTAL) 4394 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 4395 #endif 4396 #if defined(PETSC_HAVE_SCALAPACK) 4397 PETSC_INTERN PetscErrorCode MatConvert_AIJ_ScaLAPACK(Mat,MatType,MatReuse,Mat*); 4398 #endif 4399 #if defined(PETSC_HAVE_HYPRE) 4400 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 4401 #endif 4402 4403 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 4404 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 4405 PETSC_INTERN PetscErrorCode MatProductSetFromOptions_IS_XAIJ(Mat); 4406 4407 /*@C 4408 MatSeqAIJGetArray - gives read/write access to the array where the data for a MATSEQAIJ matrix is stored 4409 4410 Not Collective 4411 4412 Input Parameter: 4413 . mat - a MATSEQAIJ matrix 4414 4415 Output Parameter: 4416 . array - pointer to the data 4417 4418 Level: intermediate 4419 4420 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4421 @*/ 4422 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 4423 { 4424 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4425 4426 PetscFunctionBegin; 4427 if (aij->ops->getarray) { 4428 PetscCall((*aij->ops->getarray)(A,array)); 4429 } else { 4430 *array = aij->a; 4431 } 4432 PetscFunctionReturn(0); 4433 } 4434 4435 /*@C 4436 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 4437 4438 Not Collective 4439 4440 Input Parameters: 4441 + mat - a MATSEQAIJ matrix 4442 - array - pointer to the data 4443 4444 Level: intermediate 4445 4446 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 4447 @*/ 4448 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 4449 { 4450 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4451 4452 PetscFunctionBegin; 4453 if (aij->ops->restorearray) { 4454 PetscCall((*aij->ops->restorearray)(A,array)); 4455 } else { 4456 *array = NULL; 4457 } 4458 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4459 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4460 PetscFunctionReturn(0); 4461 } 4462 4463 /*@C 4464 MatSeqAIJGetArrayRead - gives read-only access to the array where the data for a MATSEQAIJ matrix is stored 4465 4466 Not Collective 4467 4468 Input Parameter: 4469 . mat - a MATSEQAIJ matrix 4470 4471 Output Parameter: 4472 . array - pointer to the data 4473 4474 Level: intermediate 4475 4476 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead() 4477 @*/ 4478 PetscErrorCode MatSeqAIJGetArrayRead(Mat A,const PetscScalar **array) 4479 { 4480 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4481 4482 PetscFunctionBegin; 4483 if (aij->ops->getarrayread) { 4484 PetscCall((*aij->ops->getarrayread)(A,array)); 4485 } else { 4486 *array = aij->a; 4487 } 4488 PetscFunctionReturn(0); 4489 } 4490 4491 /*@C 4492 MatSeqAIJRestoreArrayRead - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4493 4494 Not Collective 4495 4496 Input Parameter: 4497 . mat - a MATSEQAIJ matrix 4498 4499 Output Parameter: 4500 . array - pointer to the data 4501 4502 Level: intermediate 4503 4504 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4505 @*/ 4506 PetscErrorCode MatSeqAIJRestoreArrayRead(Mat A,const PetscScalar **array) 4507 { 4508 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4509 4510 PetscFunctionBegin; 4511 if (aij->ops->restorearrayread) { 4512 PetscCall((*aij->ops->restorearrayread)(A,array)); 4513 } else { 4514 *array = NULL; 4515 } 4516 PetscFunctionReturn(0); 4517 } 4518 4519 /*@C 4520 MatSeqAIJGetArrayWrite - gives write-only access to the array where the data for a MATSEQAIJ matrix is stored 4521 4522 Not Collective 4523 4524 Input Parameter: 4525 . mat - a MATSEQAIJ matrix 4526 4527 Output Parameter: 4528 . array - pointer to the data 4529 4530 Level: intermediate 4531 4532 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayRead() 4533 @*/ 4534 PetscErrorCode MatSeqAIJGetArrayWrite(Mat A,PetscScalar **array) 4535 { 4536 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4537 4538 PetscFunctionBegin; 4539 if (aij->ops->getarraywrite) { 4540 PetscCall((*aij->ops->getarraywrite)(A,array)); 4541 } else { 4542 *array = aij->a; 4543 } 4544 PetscCall(MatSeqAIJInvalidateDiagonal(A)); 4545 PetscCall(PetscObjectStateIncrease((PetscObject)A)); 4546 PetscFunctionReturn(0); 4547 } 4548 4549 /*@C 4550 MatSeqAIJRestoreArrayWrite - restore the read-only access array obtained from MatSeqAIJGetArrayRead 4551 4552 Not Collective 4553 4554 Input Parameter: 4555 . mat - a MATSEQAIJ matrix 4556 4557 Output Parameter: 4558 . array - pointer to the data 4559 4560 Level: intermediate 4561 4562 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4563 @*/ 4564 PetscErrorCode MatSeqAIJRestoreArrayWrite(Mat A,PetscScalar **array) 4565 { 4566 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4567 4568 PetscFunctionBegin; 4569 if (aij->ops->restorearraywrite) { 4570 PetscCall((*aij->ops->restorearraywrite)(A,array)); 4571 } else { 4572 *array = NULL; 4573 } 4574 PetscFunctionReturn(0); 4575 } 4576 4577 /*@C 4578 MatSeqAIJGetCSRAndMemType - Get the CSR arrays and the memory type of the SEQAIJ matrix 4579 4580 Not Collective 4581 4582 Input Parameter: 4583 . mat - a matrix of type MATSEQAIJ or its subclasses 4584 4585 Output Parameters: 4586 + i - row map array of the matrix 4587 . j - column index array of the matrix 4588 . a - data array of the matrix 4589 - memtype - memory type of the arrays 4590 4591 Notes: 4592 Any of the output parameters can be NULL, in which case the corresponding value is not returned. 4593 If mat is a device matrix, the arrays are on the device. Otherwise, they are on the host. 4594 4595 One can call this routine on a preallocated but not assembled matrix to just get the memory of the CSR underneath the matrix. 4596 If the matrix is assembled, the data array 'a' is guaranteed to have the latest values of the matrix. 4597 4598 Level: Developer 4599 4600 .seealso: MatSeqAIJGetArray(), MatSeqAIJGetArrayRead() 4601 @*/ 4602 PetscErrorCode MatSeqAIJGetCSRAndMemType(Mat mat,const PetscInt **i,const PetscInt **j,PetscScalar **a,PetscMemType *mtype) 4603 { 4604 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 4605 4606 PetscFunctionBegin; 4607 PetscCheck(mat->preallocated,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"matrix is not preallocated"); 4608 if (aij->ops->getcsrandmemtype) { 4609 PetscCall((*aij->ops->getcsrandmemtype)(mat,i,j,a,mtype)); 4610 } else { 4611 if (i) *i = aij->i; 4612 if (j) *j = aij->j; 4613 if (a) *a = aij->a; 4614 if (mtype) *mtype = PETSC_MEMTYPE_HOST; 4615 } 4616 PetscFunctionReturn(0); 4617 } 4618 4619 /*@C 4620 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 4621 4622 Not Collective 4623 4624 Input Parameter: 4625 . mat - a MATSEQAIJ matrix 4626 4627 Output Parameter: 4628 . nz - the maximum number of nonzeros in any row 4629 4630 Level: intermediate 4631 4632 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 4633 @*/ 4634 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 4635 { 4636 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 4637 4638 PetscFunctionBegin; 4639 *nz = aij->rmax; 4640 PetscFunctionReturn(0); 4641 } 4642 4643 PetscErrorCode MatSetPreallocationCOO_SeqAIJ(Mat mat, PetscCount coo_n, const PetscInt coo_i[], const PetscInt coo_j[]) 4644 { 4645 MPI_Comm comm; 4646 PetscInt *i,*j; 4647 PetscInt M,N,row; 4648 PetscCount k,p,q,nneg,nnz,start,end; /* Index the coo array, so use PetscCount as their type */ 4649 PetscInt *Ai; /* Change to PetscCount once we use it for row pointers */ 4650 PetscInt *Aj; 4651 PetscScalar *Aa; 4652 Mat_SeqAIJ *seqaij = (Mat_SeqAIJ*)(mat->data); 4653 MatType rtype; 4654 PetscCount *perm,*jmap; 4655 4656 PetscFunctionBegin; 4657 PetscCall(MatResetPreallocationCOO_SeqAIJ(mat)); 4658 PetscCall(PetscObjectGetComm((PetscObject)mat,&comm)); 4659 PetscCall(MatGetSize(mat,&M,&N)); 4660 PetscCall(PetscMalloc2(coo_n,&i,coo_n,&j)); 4661 PetscCall(PetscArraycpy(i,coo_i,coo_n)); /* Make a copy since we'll modify it */ 4662 PetscCall(PetscArraycpy(j,coo_j,coo_n)); 4663 PetscCall(PetscMalloc1(coo_n,&perm)); 4664 for (k=0; k<coo_n; k++) { /* Ignore entries with negative row or col indices */ 4665 if (j[k] < 0) i[k] = -1; 4666 perm[k] = k; 4667 } 4668 4669 /* Sort by row */ 4670 PetscCall(PetscSortIntWithIntCountArrayPair(coo_n,i,j,perm)); 4671 for (k=0; k<coo_n; k++) {if (i[k] >= 0) break;} /* Advance k to the first row with a non-negative index */ 4672 nneg = k; 4673 PetscCall(PetscMalloc1(coo_n-nneg+1,&jmap)); /* +1 to make a CSR-like data structure. jmap[i] originally is the number of repeats for i-th nonzero */ 4674 nnz = 0; /* Total number of unique nonzeros to be counted */ 4675 jmap++; /* Inc jmap by 1 for convinience */ 4676 4677 PetscCall(PetscCalloc1(M+1,&Ai)); /* CSR of A */ 4678 PetscCall(PetscMalloc1(coo_n-nneg,&Aj)); /* We have at most coo_n-nneg unique nonzeros */ 4679 4680 /* In each row, sort by column, then unique column indices to get row length */ 4681 Ai++; /* Inc by 1 for convinience */ 4682 q = 0; /* q-th unique nonzero, with q starting from 0 */ 4683 while (k<coo_n) { 4684 row = i[k]; 4685 start = k; /* [start,end) indices for this row */ 4686 while (k<coo_n && i[k] == row) k++; 4687 end = k; 4688 PetscCall(PetscSortIntWithCountArray(end-start,j+start,perm+start)); 4689 /* Find number of unique col entries in this row */ 4690 Aj[q] = j[start]; /* Log the first nonzero in this row */ 4691 jmap[q] = 1; /* Number of repeats of this nozero entry */ 4692 Ai[row] = 1; 4693 nnz++; 4694 4695 for (p=start+1; p<end; p++) { /* Scan remaining nonzero in this row */ 4696 if (j[p] != j[p-1]) { /* Meet a new nonzero */ 4697 q++; 4698 jmap[q] = 1; 4699 Aj[q] = j[p]; 4700 Ai[row]++; 4701 nnz++; 4702 } else { 4703 jmap[q]++; 4704 } 4705 } 4706 q++; /* Move to next row and thus next unique nonzero */ 4707 } 4708 PetscCall(PetscFree2(i,j)); 4709 4710 Ai--; /* Back to the beginning of Ai[] */ 4711 for (k=0; k<M; k++) Ai[k+1] += Ai[k]; 4712 jmap--; /* Back to the beginning of jmap[] */ 4713 jmap[0] = 0; 4714 for (k=0; k<nnz; k++) jmap[k+1] += jmap[k]; 4715 if (nnz < coo_n-nneg) { /* Realloc with actual number of unique nonzeros */ 4716 PetscCount *jmap_new; 4717 PetscInt *Aj_new; 4718 4719 PetscCall(PetscMalloc1(nnz+1,&jmap_new)); 4720 PetscCall(PetscArraycpy(jmap_new,jmap,nnz+1)); 4721 PetscCall(PetscFree(jmap)); 4722 jmap = jmap_new; 4723 4724 PetscCall(PetscMalloc1(nnz,&Aj_new)); 4725 PetscCall(PetscArraycpy(Aj_new,Aj,nnz)); 4726 PetscCall(PetscFree(Aj)); 4727 Aj = Aj_new; 4728 } 4729 4730 if (nneg) { /* Discard heading entries with negative indices in perm[], as we'll access it from index 0 in MatSetValuesCOO */ 4731 PetscCount *perm_new; 4732 4733 PetscCall(PetscMalloc1(coo_n-nneg,&perm_new)); 4734 PetscCall(PetscArraycpy(perm_new,perm+nneg,coo_n-nneg)); 4735 PetscCall(PetscFree(perm)); 4736 perm = perm_new; 4737 } 4738 4739 PetscCall(MatGetRootType_Private(mat,&rtype)); 4740 PetscCall(PetscCalloc1(nnz,&Aa)); /* Zero the matrix */ 4741 PetscCall(MatSetSeqAIJWithArrays_private(PETSC_COMM_SELF,M,N,Ai,Aj,Aa,rtype,mat)); 4742 4743 seqaij->singlemalloc = PETSC_FALSE; /* Ai, Aj and Aa are not allocated in one big malloc */ 4744 seqaij->free_a = seqaij->free_ij = PETSC_TRUE; /* Let newmat own Ai, Aj and Aa */ 4745 /* Record COO fields */ 4746 seqaij->coo_n = coo_n; 4747 seqaij->Atot = coo_n-nneg; /* Annz is seqaij->nz, so no need to record that again */ 4748 seqaij->jmap = jmap; /* of length nnz+1 */ 4749 seqaij->perm = perm; 4750 PetscFunctionReturn(0); 4751 } 4752 4753 static PetscErrorCode MatSetValuesCOO_SeqAIJ(Mat A,const PetscScalar v[],InsertMode imode) 4754 { 4755 Mat_SeqAIJ *aseq = (Mat_SeqAIJ*)A->data; 4756 PetscCount i,j,Annz = aseq->nz; 4757 PetscCount *perm = aseq->perm,*jmap = aseq->jmap; 4758 PetscScalar *Aa; 4759 4760 PetscFunctionBegin; 4761 PetscCall(MatSeqAIJGetArray(A,&Aa)); 4762 for (i=0; i<Annz; i++) { 4763 PetscScalar sum = 0.0; 4764 for (j=jmap[i]; j<jmap[i+1]; j++) sum += v[perm[j]]; 4765 Aa[i] = (imode == INSERT_VALUES? 0.0 : Aa[i]) + sum; 4766 } 4767 PetscCall(MatSeqAIJRestoreArray(A,&Aa)); 4768 PetscFunctionReturn(0); 4769 } 4770 4771 #if defined(PETSC_HAVE_CUDA) 4772 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCUSPARSE(Mat,MatType,MatReuse,Mat*); 4773 #endif 4774 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4775 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJKokkos(Mat,MatType,MatReuse,Mat*); 4776 #endif 4777 4778 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 4779 { 4780 Mat_SeqAIJ *b; 4781 PetscMPIInt size; 4782 4783 PetscFunctionBegin; 4784 PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)B),&size)); 4785 PetscCheckFalse(size > 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4786 4787 PetscCall(PetscNewLog(B,&b)); 4788 4789 B->data = (void*)b; 4790 4791 PetscCall(PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps))); 4792 if (B->sortedfull) B->ops->setvalues = MatSetValues_SeqAIJ_SortedFull; 4793 4794 b->row = NULL; 4795 b->col = NULL; 4796 b->icol = NULL; 4797 b->reallocs = 0; 4798 b->ignorezeroentries = PETSC_FALSE; 4799 b->roworiented = PETSC_TRUE; 4800 b->nonew = 0; 4801 b->diag = NULL; 4802 b->solve_work = NULL; 4803 B->spptr = NULL; 4804 b->saved_values = NULL; 4805 b->idiag = NULL; 4806 b->mdiag = NULL; 4807 b->ssor_work = NULL; 4808 b->omega = 1.0; 4809 b->fshift = 0.0; 4810 b->idiagvalid = PETSC_FALSE; 4811 b->ibdiagvalid = PETSC_FALSE; 4812 b->keepnonzeropattern = PETSC_FALSE; 4813 4814 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4815 4816 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4817 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ)); 4818 PetscCall(PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ)); 4819 #endif 4820 4821 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ)); 4822 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ)); 4823 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ)); 4824 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ)); 4825 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ)); 4826 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM)); 4827 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijsell_C",MatConvert_SeqAIJ_SeqAIJSELL)); 4828 #if defined(PETSC_HAVE_MKL_SPARSE) 4829 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL)); 4830 #endif 4831 #if defined(PETSC_HAVE_CUDA) 4832 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcusparse_C",MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 4833 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaijcusparse_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4834 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaijcusparse_C",MatProductSetFromOptions_SeqAIJ)); 4835 #endif 4836 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 4837 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijkokkos_C",MatConvert_SeqAIJ_SeqAIJKokkos)); 4838 #endif 4839 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL)); 4840 #if defined(PETSC_HAVE_ELEMENTAL) 4841 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental)); 4842 #endif 4843 #if defined(PETSC_HAVE_SCALAPACK) 4844 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_scalapack_C",MatConvert_AIJ_ScaLAPACK)); 4845 #endif 4846 #if defined(PETSC_HAVE_HYPRE) 4847 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE)); 4848 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_transpose_seqaij_seqaij_C",MatProductSetFromOptions_Transpose_AIJ_AIJ)); 4849 #endif 4850 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense)); 4851 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL)); 4852 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS)); 4853 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ)); 4854 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ)); 4855 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ)); 4856 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ)); 4857 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ)); 4858 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ)); 4859 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_is_seqaij_C",MatProductSetFromOptions_IS_XAIJ)); 4860 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqdense_seqaij_C",MatProductSetFromOptions_SeqDense_SeqAIJ)); 4861 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatProductSetFromOptions_seqaij_seqaij_C",MatProductSetFromOptions_SeqAIJ)); 4862 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJKron_C",MatSeqAIJKron_SeqAIJ)); 4863 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetPreallocationCOO_C",MatSetPreallocationCOO_SeqAIJ)); 4864 PetscCall(PetscObjectComposeFunction((PetscObject)B,"MatSetValuesCOO_C",MatSetValuesCOO_SeqAIJ)); 4865 PetscCall(MatCreate_SeqAIJ_Inode(B)); 4866 PetscCall(PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ)); 4867 PetscCall(MatSeqAIJSetTypeFromOptions(B)); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4868 PetscFunctionReturn(0); 4869 } 4870 4871 /* 4872 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4873 */ 4874 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4875 { 4876 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data,*a = (Mat_SeqAIJ*)A->data; 4877 PetscInt m = A->rmap->n,i; 4878 4879 PetscFunctionBegin; 4880 PetscCheckFalse(!A->assembled && cpvalues!=MAT_DO_NOT_COPY_VALUES,PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Cannot duplicate unassembled matrix"); 4881 4882 C->factortype = A->factortype; 4883 c->row = NULL; 4884 c->col = NULL; 4885 c->icol = NULL; 4886 c->reallocs = 0; 4887 4888 C->assembled = A->assembled; 4889 C->preallocated = A->preallocated; 4890 4891 if (A->preallocated) { 4892 PetscCall(PetscLayoutReference(A->rmap,&C->rmap)); 4893 PetscCall(PetscLayoutReference(A->cmap,&C->cmap)); 4894 4895 PetscCall(PetscMalloc1(m,&c->imax)); 4896 PetscCall(PetscMemcpy(c->imax,a->imax,m*sizeof(PetscInt))); 4897 PetscCall(PetscMalloc1(m,&c->ilen)); 4898 PetscCall(PetscMemcpy(c->ilen,a->ilen,m*sizeof(PetscInt))); 4899 PetscCall(PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt))); 4900 4901 /* allocate the matrix space */ 4902 if (mallocmatspace) { 4903 PetscCall(PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i)); 4904 PetscCall(PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt))); 4905 4906 c->singlemalloc = PETSC_TRUE; 4907 4908 PetscCall(PetscArraycpy(c->i,a->i,m+1)); 4909 if (m > 0) { 4910 PetscCall(PetscArraycpy(c->j,a->j,a->i[m])); 4911 if (cpvalues == MAT_COPY_VALUES) { 4912 const PetscScalar *aa; 4913 4914 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4915 PetscCall(PetscArraycpy(c->a,aa,a->i[m])); 4916 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 4917 } else { 4918 PetscCall(PetscArrayzero(c->a,a->i[m])); 4919 } 4920 } 4921 } 4922 4923 c->ignorezeroentries = a->ignorezeroentries; 4924 c->roworiented = a->roworiented; 4925 c->nonew = a->nonew; 4926 if (a->diag) { 4927 PetscCall(PetscMalloc1(m+1,&c->diag)); 4928 PetscCall(PetscMemcpy(c->diag,a->diag,m*sizeof(PetscInt))); 4929 PetscCall(PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt))); 4930 } else c->diag = NULL; 4931 4932 c->solve_work = NULL; 4933 c->saved_values = NULL; 4934 c->idiag = NULL; 4935 c->ssor_work = NULL; 4936 c->keepnonzeropattern = a->keepnonzeropattern; 4937 c->free_a = PETSC_TRUE; 4938 c->free_ij = PETSC_TRUE; 4939 4940 c->rmax = a->rmax; 4941 c->nz = a->nz; 4942 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4943 4944 c->compressedrow.use = a->compressedrow.use; 4945 c->compressedrow.nrows = a->compressedrow.nrows; 4946 if (a->compressedrow.use) { 4947 i = a->compressedrow.nrows; 4948 PetscCall(PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex)); 4949 PetscCall(PetscArraycpy(c->compressedrow.i,a->compressedrow.i,i+1)); 4950 PetscCall(PetscArraycpy(c->compressedrow.rindex,a->compressedrow.rindex,i)); 4951 } else { 4952 c->compressedrow.use = PETSC_FALSE; 4953 c->compressedrow.i = NULL; 4954 c->compressedrow.rindex = NULL; 4955 } 4956 c->nonzerorowcnt = a->nonzerorowcnt; 4957 C->nonzerostate = A->nonzerostate; 4958 4959 PetscCall(MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C)); 4960 } 4961 PetscCall(PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist)); 4962 PetscFunctionReturn(0); 4963 } 4964 4965 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4966 { 4967 PetscFunctionBegin; 4968 PetscCall(MatCreate(PetscObjectComm((PetscObject)A),B)); 4969 PetscCall(MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n)); 4970 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4971 PetscCall(MatSetBlockSizesFromMats(*B,A,A)); 4972 } 4973 PetscCall(MatSetType(*B,((PetscObject)A)->type_name)); 4974 PetscCall(MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE)); 4975 PetscFunctionReturn(0); 4976 } 4977 4978 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4979 { 4980 PetscBool isbinary, ishdf5; 4981 4982 PetscFunctionBegin; 4983 PetscValidHeaderSpecific(newMat,MAT_CLASSID,1); 4984 PetscValidHeaderSpecific(viewer,PETSC_VIEWER_CLASSID,2); 4985 /* force binary viewer to load .info file if it has not yet done so */ 4986 PetscCall(PetscViewerSetUp(viewer)); 4987 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary)); 4988 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERHDF5, &ishdf5)); 4989 if (isbinary) { 4990 PetscCall(MatLoad_SeqAIJ_Binary(newMat,viewer)); 4991 } else if (ishdf5) { 4992 #if defined(PETSC_HAVE_HDF5) 4993 PetscCall(MatLoad_AIJ_HDF5(newMat,viewer)); 4994 #else 4995 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"HDF5 not supported in this build.\nPlease reconfigure using --download-hdf5"); 4996 #endif 4997 } else { 4998 SETERRQ(PetscObjectComm((PetscObject)newMat),PETSC_ERR_SUP,"Viewer type %s not yet supported for reading %s matrices",((PetscObject)viewer)->type_name,((PetscObject)newMat)->type_name); 4999 } 5000 PetscFunctionReturn(0); 5001 } 5002 5003 PetscErrorCode MatLoad_SeqAIJ_Binary(Mat mat, PetscViewer viewer) 5004 { 5005 Mat_SeqAIJ *a = (Mat_SeqAIJ*)mat->data; 5006 PetscInt header[4],*rowlens,M,N,nz,sum,rows,cols,i; 5007 5008 PetscFunctionBegin; 5009 PetscCall(PetscViewerSetUp(viewer)); 5010 5011 /* read in matrix header */ 5012 PetscCall(PetscViewerBinaryRead(viewer,header,4,NULL,PETSC_INT)); 5013 PetscCheckFalse(header[0] != MAT_FILE_CLASSID,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Not a matrix object in file"); 5014 M = header[1]; N = header[2]; nz = header[3]; 5015 PetscCheckFalse(M < 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix row size (%" PetscInt_FMT ") in file is negative",M); 5016 PetscCheckFalse(N < 0,PetscObjectComm((PetscObject)viewer),PETSC_ERR_FILE_UNEXPECTED,"Matrix column size (%" PetscInt_FMT ") in file is negative",N); 5017 PetscCheckFalse(nz < 0,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk, cannot load as SeqAIJ"); 5018 5019 /* set block sizes from the viewer's .info file */ 5020 PetscCall(MatLoad_Binary_BlockSizes(mat,viewer)); 5021 /* set local and global sizes if not set already */ 5022 if (mat->rmap->n < 0) mat->rmap->n = M; 5023 if (mat->cmap->n < 0) mat->cmap->n = N; 5024 if (mat->rmap->N < 0) mat->rmap->N = M; 5025 if (mat->cmap->N < 0) mat->cmap->N = N; 5026 PetscCall(PetscLayoutSetUp(mat->rmap)); 5027 PetscCall(PetscLayoutSetUp(mat->cmap)); 5028 5029 /* check if the matrix sizes are correct */ 5030 PetscCall(MatGetSize(mat,&rows,&cols)); 5031 PetscCheckFalse(M != rows || N != cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different sizes (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5032 5033 /* read in row lengths */ 5034 PetscCall(PetscMalloc1(M,&rowlens)); 5035 PetscCall(PetscViewerBinaryRead(viewer,rowlens,M,NULL,PETSC_INT)); 5036 /* check if sum(rowlens) is same as nz */ 5037 sum = 0; for (i=0; i<M; i++) sum += rowlens[i]; 5038 PetscCheckFalse(sum != nz,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Inconsistent matrix data in file: nonzeros = %" PetscInt_FMT ", sum-row-lengths = %" PetscInt_FMT,nz,sum); 5039 /* preallocate and check sizes */ 5040 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(mat,0,rowlens)); 5041 PetscCall(MatGetSize(mat,&rows,&cols)); 5042 PetscCheckFalse(M != rows || N != cols,PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%" PetscInt_FMT ", %" PetscInt_FMT ") than the input matrix (%" PetscInt_FMT ", %" PetscInt_FMT ")",M,N,rows,cols); 5043 /* store row lengths */ 5044 PetscCall(PetscArraycpy(a->ilen,rowlens,M)); 5045 PetscCall(PetscFree(rowlens)); 5046 5047 /* fill in "i" row pointers */ 5048 a->i[0] = 0; for (i=0; i<M; i++) a->i[i+1] = a->i[i] + a->ilen[i]; 5049 /* read in "j" column indices */ 5050 PetscCall(PetscViewerBinaryRead(viewer,a->j,nz,NULL,PETSC_INT)); 5051 /* read in "a" nonzero values */ 5052 PetscCall(PetscViewerBinaryRead(viewer,a->a,nz,NULL,PETSC_SCALAR)); 5053 5054 PetscCall(MatAssemblyBegin(mat,MAT_FINAL_ASSEMBLY)); 5055 PetscCall(MatAssemblyEnd(mat,MAT_FINAL_ASSEMBLY)); 5056 PetscFunctionReturn(0); 5057 } 5058 5059 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 5060 { 5061 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 5062 const PetscScalar *aa,*ba; 5063 #if defined(PETSC_USE_COMPLEX) 5064 PetscInt k; 5065 #endif 5066 5067 PetscFunctionBegin; 5068 /* If the matrix dimensions are not equal,or no of nonzeros */ 5069 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 5070 *flg = PETSC_FALSE; 5071 PetscFunctionReturn(0); 5072 } 5073 5074 /* if the a->i are the same */ 5075 PetscCall(PetscArraycmp(a->i,b->i,A->rmap->n+1,flg)); 5076 if (!*flg) PetscFunctionReturn(0); 5077 5078 /* if a->j are the same */ 5079 PetscCall(PetscArraycmp(a->j,b->j,a->nz,flg)); 5080 if (!*flg) PetscFunctionReturn(0); 5081 5082 PetscCall(MatSeqAIJGetArrayRead(A,&aa)); 5083 PetscCall(MatSeqAIJGetArrayRead(B,&ba)); 5084 /* if a->a are the same */ 5085 #if defined(PETSC_USE_COMPLEX) 5086 for (k=0; k<a->nz; k++) { 5087 if (PetscRealPart(aa[k]) != PetscRealPart(ba[k]) || PetscImaginaryPart(aa[k]) != PetscImaginaryPart(ba[k])) { 5088 *flg = PETSC_FALSE; 5089 PetscFunctionReturn(0); 5090 } 5091 } 5092 #else 5093 PetscCall(PetscArraycmp(aa,ba,a->nz,flg)); 5094 #endif 5095 PetscCall(MatSeqAIJRestoreArrayRead(A,&aa)); 5096 PetscCall(MatSeqAIJRestoreArrayRead(B,&ba)); 5097 PetscFunctionReturn(0); 5098 } 5099 5100 /*@ 5101 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 5102 provided by the user. 5103 5104 Collective 5105 5106 Input Parameters: 5107 + comm - must be an MPI communicator of size 1 5108 . m - number of rows 5109 . n - number of columns 5110 . i - row indices; that is i[0] = 0, i[row] = i[row-1] + number of elements in that row of the matrix 5111 . j - column indices 5112 - a - matrix values 5113 5114 Output Parameter: 5115 . mat - the matrix 5116 5117 Level: intermediate 5118 5119 Notes: 5120 The i, j, and a arrays are not copied by this routine, the user must free these arrays 5121 once the matrix is destroyed and not before 5122 5123 You cannot set new nonzero locations into this matrix, that will generate an error. 5124 5125 The i and j indices are 0 based 5126 5127 The format which is used for the sparse matrix input, is equivalent to a 5128 row-major ordering.. i.e for the following matrix, the input data expected is 5129 as shown 5130 5131 $ 1 0 0 5132 $ 2 0 3 5133 $ 4 5 6 5134 $ 5135 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 5136 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 5137 $ v = {1,2,3,4,5,6} [size = 6] 5138 5139 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 5140 5141 @*/ 5142 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 5143 { 5144 PetscInt ii; 5145 Mat_SeqAIJ *aij; 5146 PetscInt jj; 5147 5148 PetscFunctionBegin; 5149 PetscCheckFalse(m > 0 && i[0],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 5150 PetscCall(MatCreate(comm,mat)); 5151 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5152 /* PetscCall(MatSetBlockSizes(*mat,,)); */ 5153 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5154 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,NULL)); 5155 aij = (Mat_SeqAIJ*)(*mat)->data; 5156 PetscCall(PetscMalloc1(m,&aij->imax)); 5157 PetscCall(PetscMalloc1(m,&aij->ilen)); 5158 5159 aij->i = i; 5160 aij->j = j; 5161 aij->a = a; 5162 aij->singlemalloc = PETSC_FALSE; 5163 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 5164 aij->free_a = PETSC_FALSE; 5165 aij->free_ij = PETSC_FALSE; 5166 5167 for (ii=0,aij->nonzerorowcnt=0,aij->rmax=0; ii<m; ii++) { 5168 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 5169 if (PetscDefined(USE_DEBUG)) { 5170 PetscCheckFalse(i[ii+1] - i[ii] < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative row length in i (row indices) row = %" PetscInt_FMT " length = %" PetscInt_FMT,ii,i[ii+1] - i[ii]); 5171 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 5172 PetscCheckFalse(j[jj] < j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is not sorted",jj-i[ii],j[jj],ii); 5173 PetscCheckFalse(j[jj] == j[jj-1],PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column entry number %" PetscInt_FMT " (actual column %" PetscInt_FMT ") in row %" PetscInt_FMT " is identical to previous entry",jj-i[ii],j[jj],ii); 5174 } 5175 } 5176 } 5177 if (PetscDefined(USE_DEBUG)) { 5178 for (ii=0; ii<aij->i[m]; ii++) { 5179 PetscCheckFalse(j[ii] < 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5180 PetscCheckFalse(j[ii] > n - 1,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Column index to large at location = %" PetscInt_FMT " index = %" PetscInt_FMT,ii,j[ii]); 5181 } 5182 } 5183 5184 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5185 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5186 PetscFunctionReturn(0); 5187 } 5188 5189 /*@C 5190 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 5191 provided by the user. 5192 5193 Collective 5194 5195 Input Parameters: 5196 + comm - must be an MPI communicator of size 1 5197 . m - number of rows 5198 . n - number of columns 5199 . i - row indices 5200 . j - column indices 5201 . a - matrix values 5202 . nz - number of nonzeros 5203 - idx - 0 or 1 based 5204 5205 Output Parameter: 5206 . mat - the matrix 5207 5208 Level: intermediate 5209 5210 Notes: 5211 The i and j indices are 0 based. The format which is used for the sparse matrix input, is equivalent to a row-major ordering. i.e for the following matrix, 5212 the input data expected is as shown 5213 .vb 5214 1 0 0 5215 2 0 3 5216 4 5 6 5217 5218 i = {0,1,1,2,2,2} 5219 j = {0,0,2,0,1,2} 5220 v = {1,2,3,4,5,6} 5221 .ve 5222 5223 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 5224 5225 @*/ 5226 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 5227 { 5228 PetscInt ii, *nnz, one = 1,row,col; 5229 5230 PetscFunctionBegin; 5231 PetscCall(PetscCalloc1(m,&nnz)); 5232 for (ii = 0; ii < nz; ii++) { 5233 nnz[i[ii] - !!idx] += 1; 5234 } 5235 PetscCall(MatCreate(comm,mat)); 5236 PetscCall(MatSetSizes(*mat,m,n,m,n)); 5237 PetscCall(MatSetType(*mat,MATSEQAIJ)); 5238 PetscCall(MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz)); 5239 for (ii = 0; ii < nz; ii++) { 5240 if (idx) { 5241 row = i[ii] - 1; 5242 col = j[ii] - 1; 5243 } else { 5244 row = i[ii]; 5245 col = j[ii]; 5246 } 5247 PetscCall(MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES)); 5248 } 5249 PetscCall(MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY)); 5250 PetscCall(MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY)); 5251 PetscCall(PetscFree(nnz)); 5252 PetscFunctionReturn(0); 5253 } 5254 5255 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 5256 { 5257 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 5258 5259 PetscFunctionBegin; 5260 a->idiagvalid = PETSC_FALSE; 5261 a->ibdiagvalid = PETSC_FALSE; 5262 5263 PetscCall(MatSeqAIJInvalidateDiagonal_Inode(A)); 5264 PetscFunctionReturn(0); 5265 } 5266 5267 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 5268 { 5269 PetscMPIInt size; 5270 5271 PetscFunctionBegin; 5272 PetscCallMPI(MPI_Comm_size(comm,&size)); 5273 if (size == 1) { 5274 if (scall == MAT_INITIAL_MATRIX) { 5275 PetscCall(MatDuplicate(inmat,MAT_COPY_VALUES,outmat)); 5276 } else { 5277 PetscCall(MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN)); 5278 } 5279 } else { 5280 PetscCall(MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat)); 5281 } 5282 PetscFunctionReturn(0); 5283 } 5284 5285 /* 5286 Permute A into C's *local* index space using rowemb,colemb. 5287 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 5288 of [0,m), colemb is in [0,n). 5289 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 5290 */ 5291 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 5292 { 5293 /* If making this function public, change the error returned in this function away from _PLIB. */ 5294 Mat_SeqAIJ *Baij; 5295 PetscBool seqaij; 5296 PetscInt m,n,*nz,i,j,count; 5297 PetscScalar v; 5298 const PetscInt *rowindices,*colindices; 5299 5300 PetscFunctionBegin; 5301 if (!B) PetscFunctionReturn(0); 5302 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 5303 PetscCall(PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij)); 5304 PetscCheck(seqaij,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 5305 if (rowemb) { 5306 PetscCall(ISGetLocalSize(rowemb,&m)); 5307 PetscCheckFalse(m != B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Row IS of size %" PetscInt_FMT " is incompatible with matrix row size %" PetscInt_FMT,m,B->rmap->n); 5308 } else { 5309 PetscCheckFalse(C->rmap->n != B->rmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 5310 } 5311 if (colemb) { 5312 PetscCall(ISGetLocalSize(colemb,&n)); 5313 PetscCheckFalse(n != B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Diag col IS of size %" PetscInt_FMT " is incompatible with input matrix col size %" PetscInt_FMT,n,B->cmap->n); 5314 } else { 5315 PetscCheckFalse(C->cmap->n != B->cmap->n,PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 5316 } 5317 5318 Baij = (Mat_SeqAIJ*)(B->data); 5319 if (pattern == DIFFERENT_NONZERO_PATTERN) { 5320 PetscCall(PetscMalloc1(B->rmap->n,&nz)); 5321 for (i=0; i<B->rmap->n; i++) { 5322 nz[i] = Baij->i[i+1] - Baij->i[i]; 5323 } 5324 PetscCall(MatSeqAIJSetPreallocation(C,0,nz)); 5325 PetscCall(PetscFree(nz)); 5326 } 5327 if (pattern == SUBSET_NONZERO_PATTERN) { 5328 PetscCall(MatZeroEntries(C)); 5329 } 5330 count = 0; 5331 rowindices = NULL; 5332 colindices = NULL; 5333 if (rowemb) { 5334 PetscCall(ISGetIndices(rowemb,&rowindices)); 5335 } 5336 if (colemb) { 5337 PetscCall(ISGetIndices(colemb,&colindices)); 5338 } 5339 for (i=0; i<B->rmap->n; i++) { 5340 PetscInt row; 5341 row = i; 5342 if (rowindices) row = rowindices[i]; 5343 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 5344 PetscInt col; 5345 col = Baij->j[count]; 5346 if (colindices) col = colindices[col]; 5347 v = Baij->a[count]; 5348 PetscCall(MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES)); 5349 ++count; 5350 } 5351 } 5352 /* FIXME: set C's nonzerostate correctly. */ 5353 /* Assembly for C is necessary. */ 5354 C->preallocated = PETSC_TRUE; 5355 C->assembled = PETSC_TRUE; 5356 C->was_assembled = PETSC_FALSE; 5357 PetscFunctionReturn(0); 5358 } 5359 5360 PetscFunctionList MatSeqAIJList = NULL; 5361 5362 /*@C 5363 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 5364 5365 Collective on Mat 5366 5367 Input Parameters: 5368 + mat - the matrix object 5369 - matype - matrix type 5370 5371 Options Database Key: 5372 . -mat_seqai_type <method> - for example seqaijcrl 5373 5374 Level: intermediate 5375 5376 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 5377 @*/ 5378 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 5379 { 5380 PetscBool sametype; 5381 PetscErrorCode (*r)(Mat,MatType,MatReuse,Mat*); 5382 5383 PetscFunctionBegin; 5384 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 5385 PetscCall(PetscObjectTypeCompare((PetscObject)mat,matype,&sametype)); 5386 if (sametype) PetscFunctionReturn(0); 5387 5388 PetscCall(PetscFunctionListFind(MatSeqAIJList,matype,&r)); 5389 PetscCheck(r,PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 5390 PetscCall((*r)(mat,matype,MAT_INPLACE_MATRIX,&mat)); 5391 PetscFunctionReturn(0); 5392 } 5393 5394 /*@C 5395 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 5396 5397 Not Collective 5398 5399 Input Parameters: 5400 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 5401 - function - routine to convert to subtype 5402 5403 Notes: 5404 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 5405 5406 Then, your matrix can be chosen with the procedural interface at runtime via the option 5407 $ -mat_seqaij_type my_mat 5408 5409 Level: advanced 5410 5411 .seealso: MatSeqAIJRegisterAll() 5412 5413 Level: advanced 5414 @*/ 5415 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 5416 { 5417 PetscFunctionBegin; 5418 PetscCall(MatInitializePackage()); 5419 PetscCall(PetscFunctionListAdd(&MatSeqAIJList,sname,function)); 5420 PetscFunctionReturn(0); 5421 } 5422 5423 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 5424 5425 /*@C 5426 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 5427 5428 Not Collective 5429 5430 Level: advanced 5431 5432 .seealso: MatRegisterAll(), MatSeqAIJRegister() 5433 @*/ 5434 PetscErrorCode MatSeqAIJRegisterAll(void) 5435 { 5436 PetscFunctionBegin; 5437 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 5438 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 5439 5440 PetscCall(MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL)); 5441 PetscCall(MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM)); 5442 PetscCall(MatSeqAIJRegister(MATSEQAIJSELL, MatConvert_SeqAIJ_SeqAIJSELL)); 5443 #if defined(PETSC_HAVE_MKL_SPARSE) 5444 PetscCall(MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL)); 5445 #endif 5446 #if defined(PETSC_HAVE_CUDA) 5447 PetscCall(MatSeqAIJRegister(MATSEQAIJCUSPARSE, MatConvert_SeqAIJ_SeqAIJCUSPARSE)); 5448 #endif 5449 #if defined(PETSC_HAVE_KOKKOS_KERNELS) 5450 PetscCall(MatSeqAIJRegister(MATSEQAIJKOKKOS, MatConvert_SeqAIJ_SeqAIJKokkos)); 5451 #endif 5452 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 5453 PetscCall(MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL)); 5454 #endif 5455 PetscFunctionReturn(0); 5456 } 5457 5458 /* 5459 Special version for direct calls from Fortran 5460 */ 5461 #include <petsc/private/fortranimpl.h> 5462 #if defined(PETSC_HAVE_FORTRAN_CAPS) 5463 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 5464 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 5465 #define matsetvaluesseqaij_ matsetvaluesseqaij 5466 #endif 5467 5468 /* Change these macros so can be used in void function */ 5469 5470 /* Change these macros so can be used in void function */ 5471 /* Identical to PetscCallVoid, except it assigns to *_ierr */ 5472 #undef PetscCall 5473 #define PetscCall(...) do { \ 5474 PetscErrorCode ierr_msv_mpiaij = __VA_ARGS__; \ 5475 if (PetscUnlikely(ierr_msv_mpiaij)) { \ 5476 *_ierr = PetscError(PETSC_COMM_SELF,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr_msv_mpiaij,PETSC_ERROR_REPEAT," "); \ 5477 return; \ 5478 } \ 5479 } while (0) 5480 5481 #undef SETERRQ 5482 #define SETERRQ(comm,ierr,...) do { \ 5483 *_ierr = PetscError(comm,__LINE__,PETSC_FUNCTION_NAME,__FILE__,ierr,PETSC_ERROR_INITIAL,__VA_ARGS__); \ 5484 return; \ 5485 } while (0) 5486 5487 PETSC_EXTERN void matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 5488 { 5489 Mat A = *AA; 5490 PetscInt m = *mm, n = *nn; 5491 InsertMode is = *isis; 5492 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 5493 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 5494 PetscInt *imax,*ai,*ailen; 5495 PetscInt *aj,nonew = a->nonew,lastcol = -1; 5496 MatScalar *ap,value,*aa; 5497 PetscBool ignorezeroentries = a->ignorezeroentries; 5498 PetscBool roworiented = a->roworiented; 5499 5500 PetscFunctionBegin; 5501 MatCheckPreallocated(A,1); 5502 imax = a->imax; 5503 ai = a->i; 5504 ailen = a->ilen; 5505 aj = a->j; 5506 aa = a->a; 5507 5508 for (k=0; k<m; k++) { /* loop over added rows */ 5509 row = im[k]; 5510 if (row < 0) continue; 5511 PetscCheck(row < A->rmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 5512 rp = aj + ai[row]; ap = aa + ai[row]; 5513 rmax = imax[row]; nrow = ailen[row]; 5514 low = 0; 5515 high = nrow; 5516 for (l=0; l<n; l++) { /* loop over added columns */ 5517 if (in[l] < 0) continue; 5518 PetscCheck(in[l] < A->cmap->n,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 5519 col = in[l]; 5520 if (roworiented) value = v[l + k*n]; 5521 else value = v[k + l*m]; 5522 5523 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 5524 5525 if (col <= lastcol) low = 0; 5526 else high = nrow; 5527 lastcol = col; 5528 while (high-low > 5) { 5529 t = (low+high)/2; 5530 if (rp[t] > col) high = t; 5531 else low = t; 5532 } 5533 for (i=low; i<high; i++) { 5534 if (rp[i] > col) break; 5535 if (rp[i] == col) { 5536 if (is == ADD_VALUES) ap[i] += value; 5537 else ap[i] = value; 5538 goto noinsert; 5539 } 5540 } 5541 if (value == 0.0 && ignorezeroentries) goto noinsert; 5542 if (nonew == 1) goto noinsert; 5543 PetscCheck(nonew != -1,PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 5544 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 5545 N = nrow++ - 1; a->nz++; high++; 5546 /* shift up all the later entries in this row */ 5547 for (ii=N; ii>=i; ii--) { 5548 rp[ii+1] = rp[ii]; 5549 ap[ii+1] = ap[ii]; 5550 } 5551 rp[i] = col; 5552 ap[i] = value; 5553 A->nonzerostate++; 5554 noinsert:; 5555 low = i + 1; 5556 } 5557 ailen[row] = nrow; 5558 } 5559 PetscFunctionReturnVoid(); 5560 } 5561 /* Undefining these here since they were redefined from their original definition above! No 5562 * other PETSc functions should be defined past this point, as it is impossible to recover the 5563 * original definitions */ 5564 #undef PetscCall 5565 #undef SETERRQ 5566