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