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