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