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