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