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 /* 1552 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1553 */ 1554 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1555 { 1556 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1557 PetscErrorCode ierr; 1558 PetscInt i,*diag,m = A->rmap->n; 1559 MatScalar *v = a->a; 1560 PetscScalar *idiag,*mdiag; 1561 1562 PetscFunctionBegin; 1563 if (a->idiagvalid) PetscFunctionReturn(0); 1564 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1565 diag = a->diag; 1566 if (!a->idiag) { 1567 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1568 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1569 v = a->a; 1570 } 1571 mdiag = a->mdiag; 1572 idiag = a->idiag; 1573 1574 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1575 for (i=0; i<m; i++) { 1576 mdiag[i] = v[diag[i]]; 1577 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1578 if (PetscRealPart(fshift)) { 1579 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1580 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1581 A->factorerror_zeropivot_value = 0.0; 1582 A->factorerror_zeropivot_row = i; 1583 } else SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1584 } 1585 idiag[i] = 1.0/v[diag[i]]; 1586 } 1587 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1588 } else { 1589 for (i=0; i<m; i++) { 1590 mdiag[i] = v[diag[i]]; 1591 idiag[i] = omega/(fshift + v[diag[i]]); 1592 } 1593 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1594 } 1595 a->idiagvalid = PETSC_TRUE; 1596 PetscFunctionReturn(0); 1597 } 1598 1599 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1600 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1601 { 1602 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1603 PetscScalar *x,d,sum,*t,scale; 1604 const MatScalar *v,*idiag=0,*mdiag; 1605 const PetscScalar *b, *bs,*xb, *ts; 1606 PetscErrorCode ierr; 1607 PetscInt n,m = A->rmap->n,i; 1608 const PetscInt *idx,*diag; 1609 1610 PetscFunctionBegin; 1611 its = its*lits; 1612 1613 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1614 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1615 a->fshift = fshift; 1616 a->omega = omega; 1617 1618 diag = a->diag; 1619 t = a->ssor_work; 1620 idiag = a->idiag; 1621 mdiag = a->mdiag; 1622 1623 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1624 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1625 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1626 if (flag == SOR_APPLY_UPPER) { 1627 /* apply (U + D/omega) to the vector */ 1628 bs = b; 1629 for (i=0; i<m; i++) { 1630 d = fshift + mdiag[i]; 1631 n = a->i[i+1] - diag[i] - 1; 1632 idx = a->j + diag[i] + 1; 1633 v = a->a + diag[i] + 1; 1634 sum = b[i]*d/omega; 1635 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1636 x[i] = sum; 1637 } 1638 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1639 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1640 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1641 PetscFunctionReturn(0); 1642 } 1643 1644 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1645 else if (flag & SOR_EISENSTAT) { 1646 /* Let A = L + U + D; where L is lower trianglar, 1647 U is upper triangular, E = D/omega; This routine applies 1648 1649 (L + E)^{-1} A (U + E)^{-1} 1650 1651 to a vector efficiently using Eisenstat's trick. 1652 */ 1653 scale = (2.0/omega) - 1.0; 1654 1655 /* x = (E + U)^{-1} b */ 1656 for (i=m-1; i>=0; i--) { 1657 n = a->i[i+1] - diag[i] - 1; 1658 idx = a->j + diag[i] + 1; 1659 v = a->a + diag[i] + 1; 1660 sum = b[i]; 1661 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1662 x[i] = sum*idiag[i]; 1663 } 1664 1665 /* t = b - (2*E - D)x */ 1666 v = a->a; 1667 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1668 1669 /* t = (E + L)^{-1}t */ 1670 ts = t; 1671 diag = a->diag; 1672 for (i=0; i<m; i++) { 1673 n = diag[i] - a->i[i]; 1674 idx = a->j + a->i[i]; 1675 v = a->a + a->i[i]; 1676 sum = t[i]; 1677 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1678 t[i] = sum*idiag[i]; 1679 /* x = x + t */ 1680 x[i] += t[i]; 1681 } 1682 1683 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1684 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1685 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1686 PetscFunctionReturn(0); 1687 } 1688 if (flag & SOR_ZERO_INITIAL_GUESS) { 1689 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1690 for (i=0; i<m; i++) { 1691 n = diag[i] - a->i[i]; 1692 idx = a->j + a->i[i]; 1693 v = a->a + a->i[i]; 1694 sum = b[i]; 1695 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1696 t[i] = sum; 1697 x[i] = sum*idiag[i]; 1698 } 1699 xb = t; 1700 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1701 } else xb = b; 1702 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1703 for (i=m-1; i>=0; i--) { 1704 n = a->i[i+1] - diag[i] - 1; 1705 idx = a->j + diag[i] + 1; 1706 v = a->a + diag[i] + 1; 1707 sum = xb[i]; 1708 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1709 if (xb == b) { 1710 x[i] = sum*idiag[i]; 1711 } else { 1712 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1713 } 1714 } 1715 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1716 } 1717 its--; 1718 } 1719 while (its--) { 1720 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1721 for (i=0; i<m; i++) { 1722 /* lower */ 1723 n = diag[i] - a->i[i]; 1724 idx = a->j + a->i[i]; 1725 v = a->a + a->i[i]; 1726 sum = b[i]; 1727 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1728 t[i] = sum; /* save application of the lower-triangular part */ 1729 /* upper */ 1730 n = a->i[i+1] - diag[i] - 1; 1731 idx = a->j + diag[i] + 1; 1732 v = a->a + diag[i] + 1; 1733 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1734 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1735 } 1736 xb = t; 1737 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1738 } else xb = b; 1739 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1740 for (i=m-1; i>=0; i--) { 1741 sum = xb[i]; 1742 if (xb == b) { 1743 /* whole matrix (no checkpointing available) */ 1744 n = a->i[i+1] - a->i[i]; 1745 idx = a->j + a->i[i]; 1746 v = a->a + a->i[i]; 1747 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1748 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1749 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1750 n = a->i[i+1] - diag[i] - 1; 1751 idx = a->j + diag[i] + 1; 1752 v = a->a + diag[i] + 1; 1753 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1754 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1755 } 1756 } 1757 if (xb == b) { 1758 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1759 } else { 1760 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1761 } 1762 } 1763 } 1764 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1765 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1766 PetscFunctionReturn(0); 1767 } 1768 1769 1770 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1771 { 1772 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1773 1774 PetscFunctionBegin; 1775 info->block_size = 1.0; 1776 info->nz_allocated = (double)a->maxnz; 1777 info->nz_used = (double)a->nz; 1778 info->nz_unneeded = (double)(a->maxnz - a->nz); 1779 info->assemblies = (double)A->num_ass; 1780 info->mallocs = (double)A->info.mallocs; 1781 info->memory = ((PetscObject)A)->mem; 1782 if (A->factortype) { 1783 info->fill_ratio_given = A->info.fill_ratio_given; 1784 info->fill_ratio_needed = A->info.fill_ratio_needed; 1785 info->factor_mallocs = A->info.factor_mallocs; 1786 } else { 1787 info->fill_ratio_given = 0; 1788 info->fill_ratio_needed = 0; 1789 info->factor_mallocs = 0; 1790 } 1791 PetscFunctionReturn(0); 1792 } 1793 1794 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1795 { 1796 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1797 PetscInt i,m = A->rmap->n - 1; 1798 PetscErrorCode ierr; 1799 const PetscScalar *xx; 1800 PetscScalar *bb; 1801 PetscInt d = 0; 1802 1803 PetscFunctionBegin; 1804 if (x && b) { 1805 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1806 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1807 for (i=0; i<N; i++) { 1808 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1809 bb[rows[i]] = diag*xx[rows[i]]; 1810 } 1811 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1812 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1813 } 1814 1815 if (a->keepnonzeropattern) { 1816 for (i=0; i<N; i++) { 1817 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1818 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1819 } 1820 if (diag != 0.0) { 1821 for (i=0; i<N; i++) { 1822 d = rows[i]; 1823 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); 1824 } 1825 for (i=0; i<N; i++) { 1826 a->a[a->diag[rows[i]]] = diag; 1827 } 1828 } 1829 } else { 1830 if (diag != 0.0) { 1831 for (i=0; i<N; i++) { 1832 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1833 if (a->ilen[rows[i]] > 0) { 1834 a->ilen[rows[i]] = 1; 1835 a->a[a->i[rows[i]]] = diag; 1836 a->j[a->i[rows[i]]] = rows[i]; 1837 } else { /* in case row was completely empty */ 1838 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1839 } 1840 } 1841 } else { 1842 for (i=0; i<N; i++) { 1843 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1844 a->ilen[rows[i]] = 0; 1845 } 1846 } 1847 A->nonzerostate++; 1848 } 1849 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1850 PetscFunctionReturn(0); 1851 } 1852 1853 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1854 { 1855 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1856 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1857 PetscErrorCode ierr; 1858 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1859 const PetscScalar *xx; 1860 PetscScalar *bb; 1861 1862 PetscFunctionBegin; 1863 if (x && b) { 1864 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1865 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1866 vecs = PETSC_TRUE; 1867 } 1868 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1869 for (i=0; i<N; i++) { 1870 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1871 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1872 1873 zeroed[rows[i]] = PETSC_TRUE; 1874 } 1875 for (i=0; i<A->rmap->n; i++) { 1876 if (!zeroed[i]) { 1877 for (j=a->i[i]; j<a->i[i+1]; j++) { 1878 if (zeroed[a->j[j]]) { 1879 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1880 a->a[j] = 0.0; 1881 } 1882 } 1883 } else if (vecs) bb[i] = diag*xx[i]; 1884 } 1885 if (x && b) { 1886 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1887 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1888 } 1889 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1890 if (diag != 0.0) { 1891 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1892 if (missing) { 1893 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1894 else { 1895 for (i=0; i<N; i++) { 1896 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1897 } 1898 } 1899 } else { 1900 for (i=0; i<N; i++) { 1901 a->a[a->diag[rows[i]]] = diag; 1902 } 1903 } 1904 } 1905 ierr = (*A->ops->assemblyend)(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1906 PetscFunctionReturn(0); 1907 } 1908 1909 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1910 { 1911 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1912 PetscInt *itmp; 1913 1914 PetscFunctionBegin; 1915 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1916 1917 *nz = a->i[row+1] - a->i[row]; 1918 if (v) *v = a->a + a->i[row]; 1919 if (idx) { 1920 itmp = a->j + a->i[row]; 1921 if (*nz) *idx = itmp; 1922 else *idx = 0; 1923 } 1924 PetscFunctionReturn(0); 1925 } 1926 1927 /* remove this function? */ 1928 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1929 { 1930 PetscFunctionBegin; 1931 PetscFunctionReturn(0); 1932 } 1933 1934 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1935 { 1936 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1937 MatScalar *v = a->a; 1938 PetscReal sum = 0.0; 1939 PetscErrorCode ierr; 1940 PetscInt i,j; 1941 1942 PetscFunctionBegin; 1943 if (type == NORM_FROBENIUS) { 1944 #if defined(PETSC_USE_REAL___FP16) 1945 PetscBLASInt one = 1,nz = a->nz; 1946 *nrm = BLASnrm2_(&nz,v,&one); 1947 #else 1948 for (i=0; i<a->nz; i++) { 1949 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1950 } 1951 *nrm = PetscSqrtReal(sum); 1952 #endif 1953 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1954 } else if (type == NORM_1) { 1955 PetscReal *tmp; 1956 PetscInt *jj = a->j; 1957 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 1958 *nrm = 0.0; 1959 for (j=0; j<a->nz; j++) { 1960 tmp[*jj++] += PetscAbsScalar(*v); v++; 1961 } 1962 for (j=0; j<A->cmap->n; j++) { 1963 if (tmp[j] > *nrm) *nrm = tmp[j]; 1964 } 1965 ierr = PetscFree(tmp);CHKERRQ(ierr); 1966 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1967 } else if (type == NORM_INFINITY) { 1968 *nrm = 0.0; 1969 for (j=0; j<A->rmap->n; j++) { 1970 v = a->a + a->i[j]; 1971 sum = 0.0; 1972 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1973 sum += PetscAbsScalar(*v); v++; 1974 } 1975 if (sum > *nrm) *nrm = sum; 1976 } 1977 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1978 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1979 PetscFunctionReturn(0); 1980 } 1981 1982 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1983 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1984 { 1985 PetscErrorCode ierr; 1986 PetscInt i,j,anzj; 1987 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1988 PetscInt an=A->cmap->N,am=A->rmap->N; 1989 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1990 1991 PetscFunctionBegin; 1992 /* Allocate space for symbolic transpose info and work array */ 1993 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 1994 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 1995 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 1996 1997 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 1998 /* Note: offset by 1 for fast conversion into csr format. */ 1999 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 2000 /* Form ati for csr format of A^T. */ 2001 for (i=0;i<an;i++) ati[i+1] += ati[i]; 2002 2003 /* Copy ati into atfill so we have locations of the next free space in atj */ 2004 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 2005 2006 /* Walk through A row-wise and mark nonzero entries of A^T. */ 2007 for (i=0;i<am;i++) { 2008 anzj = ai[i+1] - ai[i]; 2009 for (j=0;j<anzj;j++) { 2010 atj[atfill[*aj]] = i; 2011 atfill[*aj++] += 1; 2012 } 2013 } 2014 2015 /* Clean up temporary space and complete requests. */ 2016 ierr = PetscFree(atfill);CHKERRQ(ierr); 2017 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 2018 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2019 2020 b = (Mat_SeqAIJ*)((*B)->data); 2021 b->free_a = PETSC_FALSE; 2022 b->free_ij = PETSC_TRUE; 2023 b->nonew = 0; 2024 PetscFunctionReturn(0); 2025 } 2026 2027 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 2028 { 2029 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2030 Mat C; 2031 PetscErrorCode ierr; 2032 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 2033 MatScalar *array = a->a; 2034 2035 PetscFunctionBegin; 2036 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 2037 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 2038 2039 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 2040 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2041 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 2042 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 2043 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2044 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 2045 ierr = PetscFree(col);CHKERRQ(ierr); 2046 } else { 2047 C = *B; 2048 } 2049 2050 for (i=0; i<m; i++) { 2051 len = ai[i+1]-ai[i]; 2052 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2053 array += len; 2054 aj += len; 2055 } 2056 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2057 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2058 2059 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 2060 *B = C; 2061 } else { 2062 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2063 } 2064 PetscFunctionReturn(0); 2065 } 2066 2067 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2068 { 2069 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2070 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2071 MatScalar *va,*vb; 2072 PetscErrorCode ierr; 2073 PetscInt ma,na,mb,nb, i; 2074 2075 PetscFunctionBegin; 2076 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2077 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2078 if (ma!=nb || na!=mb) { 2079 *f = PETSC_FALSE; 2080 PetscFunctionReturn(0); 2081 } 2082 aii = aij->i; bii = bij->i; 2083 adx = aij->j; bdx = bij->j; 2084 va = aij->a; vb = bij->a; 2085 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2086 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2087 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2088 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2089 2090 *f = PETSC_TRUE; 2091 for (i=0; i<ma; i++) { 2092 while (aptr[i]<aii[i+1]) { 2093 PetscInt idc,idr; 2094 PetscScalar vc,vr; 2095 /* column/row index/value */ 2096 idc = adx[aptr[i]]; 2097 idr = bdx[bptr[idc]]; 2098 vc = va[aptr[i]]; 2099 vr = vb[bptr[idc]]; 2100 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2101 *f = PETSC_FALSE; 2102 goto done; 2103 } else { 2104 aptr[i]++; 2105 if (B || i!=idc) bptr[idc]++; 2106 } 2107 } 2108 } 2109 done: 2110 ierr = PetscFree(aptr);CHKERRQ(ierr); 2111 ierr = PetscFree(bptr);CHKERRQ(ierr); 2112 PetscFunctionReturn(0); 2113 } 2114 2115 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2116 { 2117 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2118 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2119 MatScalar *va,*vb; 2120 PetscErrorCode ierr; 2121 PetscInt ma,na,mb,nb, i; 2122 2123 PetscFunctionBegin; 2124 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2125 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2126 if (ma!=nb || na!=mb) { 2127 *f = PETSC_FALSE; 2128 PetscFunctionReturn(0); 2129 } 2130 aii = aij->i; bii = bij->i; 2131 adx = aij->j; bdx = bij->j; 2132 va = aij->a; vb = bij->a; 2133 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2134 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2135 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2136 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2137 2138 *f = PETSC_TRUE; 2139 for (i=0; i<ma; i++) { 2140 while (aptr[i]<aii[i+1]) { 2141 PetscInt idc,idr; 2142 PetscScalar vc,vr; 2143 /* column/row index/value */ 2144 idc = adx[aptr[i]]; 2145 idr = bdx[bptr[idc]]; 2146 vc = va[aptr[i]]; 2147 vr = vb[bptr[idc]]; 2148 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2149 *f = PETSC_FALSE; 2150 goto done; 2151 } else { 2152 aptr[i]++; 2153 if (B || i!=idc) bptr[idc]++; 2154 } 2155 } 2156 } 2157 done: 2158 ierr = PetscFree(aptr);CHKERRQ(ierr); 2159 ierr = PetscFree(bptr);CHKERRQ(ierr); 2160 PetscFunctionReturn(0); 2161 } 2162 2163 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2164 { 2165 PetscErrorCode ierr; 2166 2167 PetscFunctionBegin; 2168 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2169 PetscFunctionReturn(0); 2170 } 2171 2172 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2173 { 2174 PetscErrorCode ierr; 2175 2176 PetscFunctionBegin; 2177 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2178 PetscFunctionReturn(0); 2179 } 2180 2181 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2182 { 2183 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2184 const PetscScalar *l,*r; 2185 PetscScalar x; 2186 MatScalar *v; 2187 PetscErrorCode ierr; 2188 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz; 2189 const PetscInt *jj; 2190 2191 PetscFunctionBegin; 2192 if (ll) { 2193 /* The local size is used so that VecMPI can be passed to this routine 2194 by MatDiagonalScale_MPIAIJ */ 2195 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2196 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2197 ierr = VecGetArrayRead(ll,&l);CHKERRQ(ierr); 2198 v = a->a; 2199 for (i=0; i<m; i++) { 2200 x = l[i]; 2201 M = a->i[i+1] - a->i[i]; 2202 for (j=0; j<M; j++) (*v++) *= x; 2203 } 2204 ierr = VecRestoreArrayRead(ll,&l);CHKERRQ(ierr); 2205 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2206 } 2207 if (rr) { 2208 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2209 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2210 ierr = VecGetArrayRead(rr,&r);CHKERRQ(ierr); 2211 v = a->a; jj = a->j; 2212 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2213 ierr = VecRestoreArrayRead(rr,&r);CHKERRQ(ierr); 2214 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2215 } 2216 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2217 PetscFunctionReturn(0); 2218 } 2219 2220 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2221 { 2222 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2223 PetscErrorCode ierr; 2224 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2225 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2226 const PetscInt *irow,*icol; 2227 PetscInt nrows,ncols; 2228 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2229 MatScalar *a_new,*mat_a; 2230 Mat C; 2231 PetscBool stride; 2232 2233 PetscFunctionBegin; 2234 2235 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2236 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2237 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2238 2239 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2240 if (stride) { 2241 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2242 } else { 2243 first = 0; 2244 step = 0; 2245 } 2246 if (stride && step == 1) { 2247 /* special case of contiguous rows */ 2248 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2249 /* loop over new rows determining lens and starting points */ 2250 for (i=0; i<nrows; i++) { 2251 kstart = ai[irow[i]]; 2252 kend = kstart + ailen[irow[i]]; 2253 starts[i] = kstart; 2254 for (k=kstart; k<kend; k++) { 2255 if (aj[k] >= first) { 2256 starts[i] = k; 2257 break; 2258 } 2259 } 2260 sum = 0; 2261 while (k < kend) { 2262 if (aj[k++] >= first+ncols) break; 2263 sum++; 2264 } 2265 lens[i] = sum; 2266 } 2267 /* create submatrix */ 2268 if (scall == MAT_REUSE_MATRIX) { 2269 PetscInt n_cols,n_rows; 2270 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2271 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2272 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2273 C = *B; 2274 } else { 2275 PetscInt rbs,cbs; 2276 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2277 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2278 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2279 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2280 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2281 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2282 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2283 } 2284 c = (Mat_SeqAIJ*)C->data; 2285 2286 /* loop over rows inserting into submatrix */ 2287 a_new = c->a; 2288 j_new = c->j; 2289 i_new = c->i; 2290 2291 for (i=0; i<nrows; i++) { 2292 ii = starts[i]; 2293 lensi = lens[i]; 2294 for (k=0; k<lensi; k++) { 2295 *j_new++ = aj[ii+k] - first; 2296 } 2297 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2298 a_new += lensi; 2299 i_new[i+1] = i_new[i] + lensi; 2300 c->ilen[i] = lensi; 2301 } 2302 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2303 } else { 2304 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2305 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2306 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2307 for (i=0; i<ncols; i++) { 2308 #if defined(PETSC_USE_DEBUG) 2309 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); 2310 #endif 2311 smap[icol[i]] = i+1; 2312 } 2313 2314 /* determine lens of each row */ 2315 for (i=0; i<nrows; i++) { 2316 kstart = ai[irow[i]]; 2317 kend = kstart + a->ilen[irow[i]]; 2318 lens[i] = 0; 2319 for (k=kstart; k<kend; k++) { 2320 if (smap[aj[k]]) { 2321 lens[i]++; 2322 } 2323 } 2324 } 2325 /* Create and fill new matrix */ 2326 if (scall == MAT_REUSE_MATRIX) { 2327 PetscBool equal; 2328 2329 c = (Mat_SeqAIJ*)((*B)->data); 2330 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2331 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2332 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2333 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2334 C = *B; 2335 } else { 2336 PetscInt rbs,cbs; 2337 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2338 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2339 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2340 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2341 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2342 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2343 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2344 } 2345 c = (Mat_SeqAIJ*)(C->data); 2346 for (i=0; i<nrows; i++) { 2347 row = irow[i]; 2348 kstart = ai[row]; 2349 kend = kstart + a->ilen[row]; 2350 mat_i = c->i[i]; 2351 mat_j = c->j + mat_i; 2352 mat_a = c->a + mat_i; 2353 mat_ilen = c->ilen + i; 2354 for (k=kstart; k<kend; k++) { 2355 if ((tcol=smap[a->j[k]])) { 2356 *mat_j++ = tcol - 1; 2357 *mat_a++ = a->a[k]; 2358 (*mat_ilen)++; 2359 2360 } 2361 } 2362 } 2363 /* Free work space */ 2364 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2365 ierr = PetscFree(smap);CHKERRQ(ierr); 2366 ierr = PetscFree(lens);CHKERRQ(ierr); 2367 /* sort */ 2368 for (i = 0; i < nrows; i++) { 2369 PetscInt ilen; 2370 2371 mat_i = c->i[i]; 2372 mat_j = c->j + mat_i; 2373 mat_a = c->a + mat_i; 2374 ilen = c->ilen[i]; 2375 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2376 } 2377 } 2378 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2379 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2380 2381 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2382 *B = C; 2383 PetscFunctionReturn(0); 2384 } 2385 2386 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2387 { 2388 PetscErrorCode ierr; 2389 Mat B; 2390 2391 PetscFunctionBegin; 2392 if (scall == MAT_INITIAL_MATRIX) { 2393 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2394 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2395 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2396 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2397 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2398 *subMat = B; 2399 } else { 2400 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2401 } 2402 PetscFunctionReturn(0); 2403 } 2404 2405 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2406 { 2407 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2408 PetscErrorCode ierr; 2409 Mat outA; 2410 PetscBool row_identity,col_identity; 2411 2412 PetscFunctionBegin; 2413 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2414 2415 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2416 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2417 2418 outA = inA; 2419 outA->factortype = MAT_FACTOR_LU; 2420 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2421 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2422 2423 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2424 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2425 2426 a->row = row; 2427 2428 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2429 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2430 2431 a->col = col; 2432 2433 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2434 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2435 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2436 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2437 2438 if (!a->solve_work) { /* this matrix may have been factored before */ 2439 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2440 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2441 } 2442 2443 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2444 if (row_identity && col_identity) { 2445 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2446 } else { 2447 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2448 } 2449 PetscFunctionReturn(0); 2450 } 2451 2452 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2453 { 2454 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2455 PetscScalar oalpha = alpha; 2456 PetscErrorCode ierr; 2457 PetscBLASInt one = 1,bnz; 2458 2459 PetscFunctionBegin; 2460 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2461 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2462 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2463 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2464 PetscFunctionReturn(0); 2465 } 2466 2467 PetscErrorCode MatDestroySubMatrix_Private(Mat_SubSppt *submatj) 2468 { 2469 PetscErrorCode ierr; 2470 PetscInt i; 2471 2472 PetscFunctionBegin; 2473 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2474 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2475 2476 for (i=0; i<submatj->nrqr; ++i) { 2477 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2478 } 2479 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2480 2481 if (submatj->rbuf1) { 2482 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2483 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2484 } 2485 2486 for (i=0; i<submatj->nrqs; ++i) { 2487 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2488 } 2489 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2490 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2491 } 2492 2493 #if defined(PETSC_USE_CTABLE) 2494 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2495 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2496 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2497 #else 2498 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2499 #endif 2500 2501 if (!submatj->allcolumns) { 2502 #if defined(PETSC_USE_CTABLE) 2503 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2504 #else 2505 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2506 #endif 2507 } 2508 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2509 2510 ierr = PetscFree(submatj);CHKERRQ(ierr); 2511 PetscFunctionReturn(0); 2512 } 2513 2514 PetscErrorCode MatDestroySubMatrix_SeqAIJ(Mat C) 2515 { 2516 PetscErrorCode ierr; 2517 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2518 Mat_SubSppt *submatj = c->submatis1; 2519 2520 PetscFunctionBegin; 2521 ierr = submatj->destroy(C);CHKERRQ(ierr); 2522 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2523 PetscFunctionReturn(0); 2524 } 2525 2526 PetscErrorCode MatDestroySubMatrices_SeqAIJ(PetscInt n,Mat *mat[]) 2527 { 2528 PetscErrorCode ierr; 2529 PetscInt i; 2530 Mat C; 2531 Mat_SeqAIJ *c; 2532 Mat_SubSppt *submatj; 2533 2534 PetscFunctionBegin; 2535 for (i=0; i<n; i++) { 2536 C = (*mat)[i]; 2537 c = (Mat_SeqAIJ*)C->data; 2538 submatj = c->submatis1; 2539 if (submatj) { 2540 if (--((PetscObject)C)->refct <= 0) { 2541 ierr = (submatj->destroy)(C);CHKERRQ(ierr); 2542 ierr = MatDestroySubMatrix_Private(submatj);CHKERRQ(ierr); 2543 ierr = PetscLayoutDestroy(&C->rmap);CHKERRQ(ierr); 2544 ierr = PetscLayoutDestroy(&C->cmap);CHKERRQ(ierr); 2545 ierr = PetscHeaderDestroy(&C);CHKERRQ(ierr); 2546 } 2547 } else { 2548 ierr = MatDestroy(&C);CHKERRQ(ierr); 2549 } 2550 } 2551 2552 ierr = PetscFree(*mat);CHKERRQ(ierr); 2553 PetscFunctionReturn(0); 2554 } 2555 2556 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2557 { 2558 PetscErrorCode ierr; 2559 PetscInt i; 2560 2561 PetscFunctionBegin; 2562 if (scall == MAT_INITIAL_MATRIX) { 2563 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2564 } 2565 2566 for (i=0; i<n; i++) { 2567 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2568 } 2569 PetscFunctionReturn(0); 2570 } 2571 2572 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2573 { 2574 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2575 PetscErrorCode ierr; 2576 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2577 const PetscInt *idx; 2578 PetscInt start,end,*ai,*aj; 2579 PetscBT table; 2580 2581 PetscFunctionBegin; 2582 m = A->rmap->n; 2583 ai = a->i; 2584 aj = a->j; 2585 2586 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2587 2588 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2589 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2590 2591 for (i=0; i<is_max; i++) { 2592 /* Initialize the two local arrays */ 2593 isz = 0; 2594 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2595 2596 /* Extract the indices, assume there can be duplicate entries */ 2597 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2598 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2599 2600 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2601 for (j=0; j<n; ++j) { 2602 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2603 } 2604 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2605 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2606 2607 k = 0; 2608 for (j=0; j<ov; j++) { /* for each overlap */ 2609 n = isz; 2610 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2611 row = nidx[k]; 2612 start = ai[row]; 2613 end = ai[row+1]; 2614 for (l = start; l<end; l++) { 2615 val = aj[l]; 2616 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2617 } 2618 } 2619 } 2620 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2621 } 2622 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2623 ierr = PetscFree(nidx);CHKERRQ(ierr); 2624 PetscFunctionReturn(0); 2625 } 2626 2627 /* -------------------------------------------------------------- */ 2628 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2629 { 2630 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2631 PetscErrorCode ierr; 2632 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2633 const PetscInt *row,*col; 2634 PetscInt *cnew,j,*lens; 2635 IS icolp,irowp; 2636 PetscInt *cwork = NULL; 2637 PetscScalar *vwork = NULL; 2638 2639 PetscFunctionBegin; 2640 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2641 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2642 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2643 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2644 2645 /* determine lengths of permuted rows */ 2646 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2647 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2648 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2649 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2650 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2651 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2652 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2653 ierr = PetscFree(lens);CHKERRQ(ierr); 2654 2655 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2656 for (i=0; i<m; i++) { 2657 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2658 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2659 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2660 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2661 } 2662 ierr = PetscFree(cnew);CHKERRQ(ierr); 2663 2664 (*B)->assembled = PETSC_FALSE; 2665 2666 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2667 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2668 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2669 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2670 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2671 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2672 PetscFunctionReturn(0); 2673 } 2674 2675 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2676 { 2677 PetscErrorCode ierr; 2678 2679 PetscFunctionBegin; 2680 /* If the two matrices have the same copy implementation, use fast copy. */ 2681 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2682 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2683 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2684 2685 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"); 2686 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2687 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2688 } else { 2689 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2690 } 2691 PetscFunctionReturn(0); 2692 } 2693 2694 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2695 { 2696 PetscErrorCode ierr; 2697 2698 PetscFunctionBegin; 2699 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2700 PetscFunctionReturn(0); 2701 } 2702 2703 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2704 { 2705 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2706 2707 PetscFunctionBegin; 2708 *array = a->a; 2709 PetscFunctionReturn(0); 2710 } 2711 2712 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2713 { 2714 PetscFunctionBegin; 2715 PetscFunctionReturn(0); 2716 } 2717 2718 /* 2719 Computes the number of nonzeros per row needed for preallocation when X and Y 2720 have different nonzero structure. 2721 */ 2722 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2723 { 2724 PetscInt i,j,k,nzx,nzy; 2725 2726 PetscFunctionBegin; 2727 /* Set the number of nonzeros in the new matrix */ 2728 for (i=0; i<m; i++) { 2729 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2730 nzx = xi[i+1] - xi[i]; 2731 nzy = yi[i+1] - yi[i]; 2732 nnz[i] = 0; 2733 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2734 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2735 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2736 nnz[i]++; 2737 } 2738 for (; k<nzy; k++) nnz[i]++; 2739 } 2740 PetscFunctionReturn(0); 2741 } 2742 2743 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2744 { 2745 PetscInt m = Y->rmap->N; 2746 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2747 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2748 PetscErrorCode ierr; 2749 2750 PetscFunctionBegin; 2751 /* Set the number of nonzeros in the new matrix */ 2752 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2753 PetscFunctionReturn(0); 2754 } 2755 2756 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2757 { 2758 PetscErrorCode ierr; 2759 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2760 PetscBLASInt one=1,bnz; 2761 2762 PetscFunctionBegin; 2763 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2764 if (str == SAME_NONZERO_PATTERN) { 2765 PetscScalar alpha = a; 2766 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2767 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2768 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2769 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2770 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2771 } else { 2772 Mat B; 2773 PetscInt *nnz; 2774 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2775 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2776 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2777 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2778 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2779 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2780 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2781 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2782 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2783 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2784 ierr = PetscFree(nnz);CHKERRQ(ierr); 2785 } 2786 PetscFunctionReturn(0); 2787 } 2788 2789 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2790 { 2791 #if defined(PETSC_USE_COMPLEX) 2792 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2793 PetscInt i,nz; 2794 PetscScalar *a; 2795 2796 PetscFunctionBegin; 2797 nz = aij->nz; 2798 a = aij->a; 2799 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2800 #else 2801 PetscFunctionBegin; 2802 #endif 2803 PetscFunctionReturn(0); 2804 } 2805 2806 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2807 { 2808 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2809 PetscErrorCode ierr; 2810 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2811 PetscReal atmp; 2812 PetscScalar *x; 2813 MatScalar *aa; 2814 2815 PetscFunctionBegin; 2816 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2817 aa = a->a; 2818 ai = a->i; 2819 aj = a->j; 2820 2821 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2822 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2823 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2824 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2825 for (i=0; i<m; i++) { 2826 ncols = ai[1] - ai[0]; ai++; 2827 x[i] = 0.0; 2828 for (j=0; j<ncols; j++) { 2829 atmp = PetscAbsScalar(*aa); 2830 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2831 aa++; aj++; 2832 } 2833 } 2834 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2835 PetscFunctionReturn(0); 2836 } 2837 2838 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2839 { 2840 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2841 PetscErrorCode ierr; 2842 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2843 PetscScalar *x; 2844 MatScalar *aa; 2845 2846 PetscFunctionBegin; 2847 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2848 aa = a->a; 2849 ai = a->i; 2850 aj = a->j; 2851 2852 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2853 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2854 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2855 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2856 for (i=0; i<m; i++) { 2857 ncols = ai[1] - ai[0]; ai++; 2858 if (ncols == A->cmap->n) { /* row is dense */ 2859 x[i] = *aa; if (idx) idx[i] = 0; 2860 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2861 x[i] = 0.0; 2862 if (idx) { 2863 idx[i] = 0; /* in case ncols is zero */ 2864 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2865 if (aj[j] > j) { 2866 idx[i] = j; 2867 break; 2868 } 2869 } 2870 } 2871 } 2872 for (j=0; j<ncols; j++) { 2873 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2874 aa++; aj++; 2875 } 2876 } 2877 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2878 PetscFunctionReturn(0); 2879 } 2880 2881 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2882 { 2883 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2884 PetscErrorCode ierr; 2885 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2886 PetscReal atmp; 2887 PetscScalar *x; 2888 MatScalar *aa; 2889 2890 PetscFunctionBegin; 2891 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2892 aa = a->a; 2893 ai = a->i; 2894 aj = a->j; 2895 2896 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2897 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2898 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2899 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); 2900 for (i=0; i<m; i++) { 2901 ncols = ai[1] - ai[0]; ai++; 2902 if (ncols) { 2903 /* Get first nonzero */ 2904 for (j = 0; j < ncols; j++) { 2905 atmp = PetscAbsScalar(aa[j]); 2906 if (atmp > 1.0e-12) { 2907 x[i] = atmp; 2908 if (idx) idx[i] = aj[j]; 2909 break; 2910 } 2911 } 2912 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2913 } else { 2914 x[i] = 0.0; if (idx) idx[i] = 0; 2915 } 2916 for (j = 0; j < ncols; j++) { 2917 atmp = PetscAbsScalar(*aa); 2918 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2919 aa++; aj++; 2920 } 2921 } 2922 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2923 PetscFunctionReturn(0); 2924 } 2925 2926 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2927 { 2928 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2929 PetscErrorCode ierr; 2930 PetscInt i,j,m = A->rmap->n,ncols,n; 2931 const PetscInt *ai,*aj; 2932 PetscScalar *x; 2933 const MatScalar *aa; 2934 2935 PetscFunctionBegin; 2936 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2937 aa = a->a; 2938 ai = a->i; 2939 aj = a->j; 2940 2941 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2942 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2943 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2944 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2945 for (i=0; i<m; i++) { 2946 ncols = ai[1] - ai[0]; ai++; 2947 if (ncols == A->cmap->n) { /* row is dense */ 2948 x[i] = *aa; if (idx) idx[i] = 0; 2949 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2950 x[i] = 0.0; 2951 if (idx) { /* find first implicit 0.0 in the row */ 2952 idx[i] = 0; /* in case ncols is zero */ 2953 for (j=0; j<ncols; j++) { 2954 if (aj[j] > j) { 2955 idx[i] = j; 2956 break; 2957 } 2958 } 2959 } 2960 } 2961 for (j=0; j<ncols; j++) { 2962 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2963 aa++; aj++; 2964 } 2965 } 2966 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2967 PetscFunctionReturn(0); 2968 } 2969 2970 #include <petscblaslapack.h> 2971 #include <petsc/private/kernels/blockinvert.h> 2972 2973 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2974 { 2975 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2976 PetscErrorCode ierr; 2977 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2978 MatScalar *diag,work[25],*v_work; 2979 PetscReal shift = 0.0; 2980 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2981 2982 PetscFunctionBegin; 2983 allowzeropivot = PetscNot(A->erroriffailure); 2984 if (a->ibdiagvalid) { 2985 if (values) *values = a->ibdiag; 2986 PetscFunctionReturn(0); 2987 } 2988 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2989 if (!a->ibdiag) { 2990 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2991 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2992 } 2993 diag = a->ibdiag; 2994 if (values) *values = a->ibdiag; 2995 /* factor and invert each block */ 2996 switch (bs) { 2997 case 1: 2998 for (i=0; i<mbs; i++) { 2999 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 3000 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 3001 if (allowzeropivot) { 3002 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3003 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 3004 A->factorerror_zeropivot_row = i; 3005 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 3006 } 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); 3007 } 3008 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 3009 } 3010 break; 3011 case 2: 3012 for (i=0; i<mbs; i++) { 3013 ij[0] = 2*i; ij[1] = 2*i + 1; 3014 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 3015 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3016 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3017 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 3018 diag += 4; 3019 } 3020 break; 3021 case 3: 3022 for (i=0; i<mbs; i++) { 3023 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 3024 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 3025 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3026 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3027 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 3028 diag += 9; 3029 } 3030 break; 3031 case 4: 3032 for (i=0; i<mbs; i++) { 3033 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 3034 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 3035 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3036 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3037 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 3038 diag += 16; 3039 } 3040 break; 3041 case 5: 3042 for (i=0; i<mbs; i++) { 3043 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 3044 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3045 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3046 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3047 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3048 diag += 25; 3049 } 3050 break; 3051 case 6: 3052 for (i=0; i<mbs; i++) { 3053 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; 3054 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3055 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3056 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3057 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3058 diag += 36; 3059 } 3060 break; 3061 case 7: 3062 for (i=0; i<mbs; i++) { 3063 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; 3064 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3065 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3066 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3067 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3068 diag += 49; 3069 } 3070 break; 3071 default: 3072 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3073 for (i=0; i<mbs; i++) { 3074 for (j=0; j<bs; j++) { 3075 IJ[j] = bs*i + j; 3076 } 3077 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3078 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3079 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3080 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3081 diag += bs2; 3082 } 3083 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3084 } 3085 a->ibdiagvalid = PETSC_TRUE; 3086 PetscFunctionReturn(0); 3087 } 3088 3089 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3090 { 3091 PetscErrorCode ierr; 3092 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3093 PetscScalar a; 3094 PetscInt m,n,i,j,col; 3095 3096 PetscFunctionBegin; 3097 if (!x->assembled) { 3098 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3099 for (i=0; i<m; i++) { 3100 for (j=0; j<aij->imax[i]; j++) { 3101 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3102 col = (PetscInt)(n*PetscRealPart(a)); 3103 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3104 } 3105 } 3106 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3107 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3108 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3109 PetscFunctionReturn(0); 3110 } 3111 3112 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3113 { 3114 PetscErrorCode ierr; 3115 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3116 3117 PetscFunctionBegin; 3118 if (!Y->preallocated || !aij->nz) { 3119 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3120 } 3121 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3122 PetscFunctionReturn(0); 3123 } 3124 3125 /* -------------------------------------------------------------------*/ 3126 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3127 MatGetRow_SeqAIJ, 3128 MatRestoreRow_SeqAIJ, 3129 MatMult_SeqAIJ, 3130 /* 4*/ MatMultAdd_SeqAIJ, 3131 MatMultTranspose_SeqAIJ, 3132 MatMultTransposeAdd_SeqAIJ, 3133 0, 3134 0, 3135 0, 3136 /* 10*/ 0, 3137 MatLUFactor_SeqAIJ, 3138 0, 3139 MatSOR_SeqAIJ, 3140 MatTranspose_SeqAIJ, 3141 /*1 5*/ MatGetInfo_SeqAIJ, 3142 MatEqual_SeqAIJ, 3143 MatGetDiagonal_SeqAIJ, 3144 MatDiagonalScale_SeqAIJ, 3145 MatNorm_SeqAIJ, 3146 /* 20*/ 0, 3147 MatAssemblyEnd_SeqAIJ, 3148 MatSetOption_SeqAIJ, 3149 MatZeroEntries_SeqAIJ, 3150 /* 24*/ MatZeroRows_SeqAIJ, 3151 0, 3152 0, 3153 0, 3154 0, 3155 /* 29*/ MatSetUp_SeqAIJ, 3156 0, 3157 0, 3158 0, 3159 0, 3160 /* 34*/ MatDuplicate_SeqAIJ, 3161 0, 3162 0, 3163 MatILUFactor_SeqAIJ, 3164 0, 3165 /* 39*/ MatAXPY_SeqAIJ, 3166 MatCreateSubMatrices_SeqAIJ, 3167 MatIncreaseOverlap_SeqAIJ, 3168 MatGetValues_SeqAIJ, 3169 MatCopy_SeqAIJ, 3170 /* 44*/ MatGetRowMax_SeqAIJ, 3171 MatScale_SeqAIJ, 3172 MatShift_SeqAIJ, 3173 MatDiagonalSet_SeqAIJ, 3174 MatZeroRowsColumns_SeqAIJ, 3175 /* 49*/ MatSetRandom_SeqAIJ, 3176 MatGetRowIJ_SeqAIJ, 3177 MatRestoreRowIJ_SeqAIJ, 3178 MatGetColumnIJ_SeqAIJ, 3179 MatRestoreColumnIJ_SeqAIJ, 3180 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3181 0, 3182 0, 3183 MatPermute_SeqAIJ, 3184 0, 3185 /* 59*/ 0, 3186 MatDestroy_SeqAIJ, 3187 MatView_SeqAIJ, 3188 0, 3189 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3190 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3191 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3192 0, 3193 0, 3194 0, 3195 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3196 MatGetRowMinAbs_SeqAIJ, 3197 0, 3198 0, 3199 0, 3200 /* 74*/ 0, 3201 MatFDColoringApply_AIJ, 3202 0, 3203 0, 3204 0, 3205 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3206 0, 3207 0, 3208 0, 3209 MatLoad_SeqAIJ, 3210 /* 84*/ MatIsSymmetric_SeqAIJ, 3211 MatIsHermitian_SeqAIJ, 3212 0, 3213 0, 3214 0, 3215 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3216 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3217 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3218 MatPtAP_SeqAIJ_SeqAIJ, 3219 MatPtAPSymbolic_SeqAIJ_SeqAIJ_SparseAxpy, 3220 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ_SparseAxpy, 3221 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3222 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3223 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3224 0, 3225 /* 99*/ 0, 3226 0, 3227 0, 3228 MatConjugate_SeqAIJ, 3229 0, 3230 /*104*/ MatSetValuesRow_SeqAIJ, 3231 MatRealPart_SeqAIJ, 3232 MatImaginaryPart_SeqAIJ, 3233 0, 3234 0, 3235 /*109*/ MatMatSolve_SeqAIJ, 3236 0, 3237 MatGetRowMin_SeqAIJ, 3238 0, 3239 MatMissingDiagonal_SeqAIJ, 3240 /*114*/ 0, 3241 0, 3242 0, 3243 0, 3244 0, 3245 /*119*/ 0, 3246 0, 3247 0, 3248 0, 3249 MatGetMultiProcBlock_SeqAIJ, 3250 /*124*/ MatFindNonzeroRows_SeqAIJ, 3251 MatGetColumnNorms_SeqAIJ, 3252 MatInvertBlockDiagonal_SeqAIJ, 3253 0, 3254 0, 3255 /*129*/ 0, 3256 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3257 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3258 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3259 MatTransposeColoringCreate_SeqAIJ, 3260 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3261 MatTransColoringApplyDenToSp_SeqAIJ, 3262 MatRARt_SeqAIJ_SeqAIJ, 3263 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3264 MatRARtNumeric_SeqAIJ_SeqAIJ, 3265 /*139*/0, 3266 0, 3267 0, 3268 MatFDColoringSetUp_SeqXAIJ, 3269 MatFindOffBlockDiagonalEntries_SeqAIJ, 3270 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ, 3271 MatDestroySubMatrices_SeqAIJ 3272 }; 3273 3274 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3275 { 3276 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3277 PetscInt i,nz,n; 3278 3279 PetscFunctionBegin; 3280 nz = aij->maxnz; 3281 n = mat->rmap->n; 3282 for (i=0; i<nz; i++) { 3283 aij->j[i] = indices[i]; 3284 } 3285 aij->nz = nz; 3286 for (i=0; i<n; i++) { 3287 aij->ilen[i] = aij->imax[i]; 3288 } 3289 PetscFunctionReturn(0); 3290 } 3291 3292 /*@ 3293 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3294 in the matrix. 3295 3296 Input Parameters: 3297 + mat - the SeqAIJ matrix 3298 - indices - the column indices 3299 3300 Level: advanced 3301 3302 Notes: 3303 This can be called if you have precomputed the nonzero structure of the 3304 matrix and want to provide it to the matrix object to improve the performance 3305 of the MatSetValues() operation. 3306 3307 You MUST have set the correct numbers of nonzeros per row in the call to 3308 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3309 3310 MUST be called before any calls to MatSetValues(); 3311 3312 The indices should start with zero, not one. 3313 3314 @*/ 3315 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3316 { 3317 PetscErrorCode ierr; 3318 3319 PetscFunctionBegin; 3320 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3321 PetscValidPointer(indices,2); 3322 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3323 PetscFunctionReturn(0); 3324 } 3325 3326 /* ----------------------------------------------------------------------------------------*/ 3327 3328 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3329 { 3330 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3331 PetscErrorCode ierr; 3332 size_t nz = aij->i[mat->rmap->n]; 3333 3334 PetscFunctionBegin; 3335 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3336 3337 /* allocate space for values if not already there */ 3338 if (!aij->saved_values) { 3339 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3340 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3341 } 3342 3343 /* copy values over */ 3344 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3345 PetscFunctionReturn(0); 3346 } 3347 3348 /*@ 3349 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3350 example, reuse of the linear part of a Jacobian, while recomputing the 3351 nonlinear portion. 3352 3353 Collect on Mat 3354 3355 Input Parameters: 3356 . mat - the matrix (currently only AIJ matrices support this option) 3357 3358 Level: advanced 3359 3360 Common Usage, with SNESSolve(): 3361 $ Create Jacobian matrix 3362 $ Set linear terms into matrix 3363 $ Apply boundary conditions to matrix, at this time matrix must have 3364 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3365 $ boundary conditions again will not change the nonzero structure 3366 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3367 $ ierr = MatStoreValues(mat); 3368 $ Call SNESSetJacobian() with matrix 3369 $ In your Jacobian routine 3370 $ ierr = MatRetrieveValues(mat); 3371 $ Set nonlinear terms in matrix 3372 3373 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3374 $ // build linear portion of Jacobian 3375 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3376 $ ierr = MatStoreValues(mat); 3377 $ loop over nonlinear iterations 3378 $ ierr = MatRetrieveValues(mat); 3379 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3380 $ // call MatAssemblyBegin/End() on matrix 3381 $ Solve linear system with Jacobian 3382 $ endloop 3383 3384 Notes: 3385 Matrix must already be assemblied before calling this routine 3386 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3387 calling this routine. 3388 3389 When this is called multiple times it overwrites the previous set of stored values 3390 and does not allocated additional space. 3391 3392 .seealso: MatRetrieveValues() 3393 3394 @*/ 3395 PetscErrorCode MatStoreValues(Mat mat) 3396 { 3397 PetscErrorCode ierr; 3398 3399 PetscFunctionBegin; 3400 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3401 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3402 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3403 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3404 PetscFunctionReturn(0); 3405 } 3406 3407 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3408 { 3409 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3410 PetscErrorCode ierr; 3411 PetscInt nz = aij->i[mat->rmap->n]; 3412 3413 PetscFunctionBegin; 3414 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3415 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3416 /* copy values over */ 3417 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3418 PetscFunctionReturn(0); 3419 } 3420 3421 /*@ 3422 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3423 example, reuse of the linear part of a Jacobian, while recomputing the 3424 nonlinear portion. 3425 3426 Collect on Mat 3427 3428 Input Parameters: 3429 . mat - the matrix (currently only AIJ matrices support this option) 3430 3431 Level: advanced 3432 3433 .seealso: MatStoreValues() 3434 3435 @*/ 3436 PetscErrorCode MatRetrieveValues(Mat mat) 3437 { 3438 PetscErrorCode ierr; 3439 3440 PetscFunctionBegin; 3441 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3442 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3443 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3444 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3445 PetscFunctionReturn(0); 3446 } 3447 3448 3449 /* --------------------------------------------------------------------------------*/ 3450 /*@C 3451 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3452 (the default parallel PETSc format). For good matrix assembly performance 3453 the user should preallocate the matrix storage by setting the parameter nz 3454 (or the array nnz). By setting these parameters accurately, performance 3455 during matrix assembly can be increased by more than a factor of 50. 3456 3457 Collective on MPI_Comm 3458 3459 Input Parameters: 3460 + comm - MPI communicator, set to PETSC_COMM_SELF 3461 . m - number of rows 3462 . n - number of columns 3463 . nz - number of nonzeros per row (same for all rows) 3464 - nnz - array containing the number of nonzeros in the various rows 3465 (possibly different for each row) or NULL 3466 3467 Output Parameter: 3468 . A - the matrix 3469 3470 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3471 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3472 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3473 3474 Notes: 3475 If nnz is given then nz is ignored 3476 3477 The AIJ format (also called the Yale sparse matrix format or 3478 compressed row storage), is fully compatible with standard Fortran 77 3479 storage. That is, the stored row and column indices can begin at 3480 either one (as in Fortran) or zero. See the users' manual for details. 3481 3482 Specify the preallocated storage with either nz or nnz (not both). 3483 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3484 allocation. For large problems you MUST preallocate memory or you 3485 will get TERRIBLE performance, see the users' manual chapter on matrices. 3486 3487 By default, this format uses inodes (identical nodes) when possible, to 3488 improve numerical efficiency of matrix-vector products and solves. We 3489 search for consecutive rows with the same nonzero structure, thereby 3490 reusing matrix information to achieve increased efficiency. 3491 3492 Options Database Keys: 3493 + -mat_no_inode - Do not use inodes 3494 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3495 3496 Level: intermediate 3497 3498 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3499 3500 @*/ 3501 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3502 { 3503 PetscErrorCode ierr; 3504 3505 PetscFunctionBegin; 3506 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3507 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3508 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3509 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3510 PetscFunctionReturn(0); 3511 } 3512 3513 /*@C 3514 MatSeqAIJSetPreallocation - For good matrix assembly performance 3515 the user should preallocate the matrix storage by setting the parameter nz 3516 (or the array nnz). By setting these parameters accurately, performance 3517 during matrix assembly can be increased by more than a factor of 50. 3518 3519 Collective on MPI_Comm 3520 3521 Input Parameters: 3522 + B - The matrix 3523 . nz - number of nonzeros per row (same for all rows) 3524 - nnz - array containing the number of nonzeros in the various rows 3525 (possibly different for each row) or NULL 3526 3527 Notes: 3528 If nnz is given then nz is ignored 3529 3530 The AIJ format (also called the Yale sparse matrix format or 3531 compressed row storage), is fully compatible with standard Fortran 77 3532 storage. That is, the stored row and column indices can begin at 3533 either one (as in Fortran) or zero. See the users' manual for details. 3534 3535 Specify the preallocated storage with either nz or nnz (not both). 3536 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3537 allocation. For large problems you MUST preallocate memory or you 3538 will get TERRIBLE performance, see the users' manual chapter on matrices. 3539 3540 You can call MatGetInfo() to get information on how effective the preallocation was; 3541 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3542 You can also run with the option -info and look for messages with the string 3543 malloc in them to see if additional memory allocation was needed. 3544 3545 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3546 entries or columns indices 3547 3548 By default, this format uses inodes (identical nodes) when possible, to 3549 improve numerical efficiency of matrix-vector products and solves. We 3550 search for consecutive rows with the same nonzero structure, thereby 3551 reusing matrix information to achieve increased efficiency. 3552 3553 Options Database Keys: 3554 + -mat_no_inode - Do not use inodes 3555 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3556 3557 Level: intermediate 3558 3559 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3560 3561 @*/ 3562 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3563 { 3564 PetscErrorCode ierr; 3565 3566 PetscFunctionBegin; 3567 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3568 PetscValidType(B,1); 3569 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3570 PetscFunctionReturn(0); 3571 } 3572 3573 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3574 { 3575 Mat_SeqAIJ *b; 3576 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3577 PetscErrorCode ierr; 3578 PetscInt i; 3579 3580 PetscFunctionBegin; 3581 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3582 if (nz == MAT_SKIP_ALLOCATION) { 3583 skipallocation = PETSC_TRUE; 3584 nz = 0; 3585 } 3586 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3587 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3588 3589 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3590 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3591 if (nnz) { 3592 for (i=0; i<B->rmap->n; i++) { 3593 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]); 3594 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); 3595 } 3596 } 3597 3598 B->preallocated = PETSC_TRUE; 3599 3600 b = (Mat_SeqAIJ*)B->data; 3601 3602 if (!skipallocation) { 3603 if (!b->imax) { 3604 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3605 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3606 } 3607 if (!b->ipre) { 3608 ierr = PetscMalloc1(B->rmap->n,&b->ipre);CHKERRQ(ierr); 3609 ierr = PetscLogObjectMemory((PetscObject)B,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3610 } 3611 if (!nnz) { 3612 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3613 else if (nz < 0) nz = 1; 3614 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3615 nz = nz*B->rmap->n; 3616 } else { 3617 nz = 0; 3618 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3619 } 3620 /* b->ilen will count nonzeros in each row so far. */ 3621 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3622 3623 /* allocate the matrix space */ 3624 /* FIXME: should B's old memory be unlogged? */ 3625 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3626 if (B->structure_only) { 3627 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3628 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3629 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3630 } else { 3631 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3632 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3633 } 3634 b->i[0] = 0; 3635 for (i=1; i<B->rmap->n+1; i++) { 3636 b->i[i] = b->i[i-1] + b->imax[i-1]; 3637 } 3638 if (B->structure_only) { 3639 b->singlemalloc = PETSC_FALSE; 3640 b->free_a = PETSC_FALSE; 3641 } else { 3642 b->singlemalloc = PETSC_TRUE; 3643 b->free_a = PETSC_TRUE; 3644 } 3645 b->free_ij = PETSC_TRUE; 3646 } else { 3647 b->free_a = PETSC_FALSE; 3648 b->free_ij = PETSC_FALSE; 3649 } 3650 3651 if (b->ipre && nnz != b->ipre && b->imax) { 3652 /* reserve user-requested sparsity */ 3653 ierr = PetscMemcpy(b->ipre,b->imax,B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3654 } 3655 3656 3657 b->nz = 0; 3658 b->maxnz = nz; 3659 B->info.nz_unneeded = (double)b->maxnz; 3660 if (realalloc) { 3661 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3662 } 3663 B->was_assembled = PETSC_FALSE; 3664 B->assembled = PETSC_FALSE; 3665 PetscFunctionReturn(0); 3666 } 3667 3668 3669 PetscErrorCode MatResetPreallocation_SeqAIJ(Mat A) 3670 { 3671 Mat_SeqAIJ *a; 3672 PetscInt i; 3673 PetscErrorCode ierr; 3674 3675 PetscFunctionBegin; 3676 PetscValidHeaderSpecific(A,MAT_CLASSID,1); 3677 a = (Mat_SeqAIJ*)A->data; 3678 /* if no saved info, we error out */ 3679 if (!a->ipre) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_ARG_NULL,"No saved preallocation info \n"); 3680 3681 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"); 3682 3683 ierr = PetscMemcpy(a->imax,a->ipre,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3684 ierr = PetscMemzero(a->ilen,A->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3685 a->i[0] = 0; 3686 for (i=1; i<A->rmap->n+1; i++) { 3687 a->i[i] = a->i[i-1] + a->imax[i-1]; 3688 } 3689 A->preallocated = PETSC_TRUE; 3690 a->nz = 0; 3691 a->maxnz = a->i[A->rmap->n]; 3692 A->info.nz_unneeded = (double)a->maxnz; 3693 A->was_assembled = PETSC_FALSE; 3694 A->assembled = PETSC_FALSE; 3695 PetscFunctionReturn(0); 3696 } 3697 3698 /*@ 3699 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3700 3701 Input Parameters: 3702 + B - the matrix 3703 . i - the indices into j for the start of each row (starts with zero) 3704 . j - the column indices for each row (starts with zero) these must be sorted for each row 3705 - v - optional values in the matrix 3706 3707 Level: developer 3708 3709 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3710 3711 .keywords: matrix, aij, compressed row, sparse, sequential 3712 3713 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), MATSEQAIJ 3714 @*/ 3715 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3716 { 3717 PetscErrorCode ierr; 3718 3719 PetscFunctionBegin; 3720 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3721 PetscValidType(B,1); 3722 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3723 PetscFunctionReturn(0); 3724 } 3725 3726 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3727 { 3728 PetscInt i; 3729 PetscInt m,n; 3730 PetscInt nz; 3731 PetscInt *nnz, nz_max = 0; 3732 PetscScalar *values; 3733 PetscErrorCode ierr; 3734 3735 PetscFunctionBegin; 3736 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3737 3738 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3739 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3740 3741 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3742 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3743 for (i = 0; i < m; i++) { 3744 nz = Ii[i+1]- Ii[i]; 3745 nz_max = PetscMax(nz_max, nz); 3746 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3747 nnz[i] = nz; 3748 } 3749 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3750 ierr = PetscFree(nnz);CHKERRQ(ierr); 3751 3752 if (v) { 3753 values = (PetscScalar*) v; 3754 } else { 3755 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3756 } 3757 3758 for (i = 0; i < m; i++) { 3759 nz = Ii[i+1] - Ii[i]; 3760 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3761 } 3762 3763 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3764 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3765 3766 if (!v) { 3767 ierr = PetscFree(values);CHKERRQ(ierr); 3768 } 3769 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3770 PetscFunctionReturn(0); 3771 } 3772 3773 #include <../src/mat/impls/dense/seq/dense.h> 3774 #include <petsc/private/kernels/petscaxpy.h> 3775 3776 /* 3777 Computes (B'*A')' since computing B*A directly is untenable 3778 3779 n p p 3780 ( ) ( ) ( ) 3781 m ( A ) * n ( B ) = m ( C ) 3782 ( ) ( ) ( ) 3783 3784 */ 3785 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3786 { 3787 PetscErrorCode ierr; 3788 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3789 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3790 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3791 PetscInt i,n,m,q,p; 3792 const PetscInt *ii,*idx; 3793 const PetscScalar *b,*a,*a_q; 3794 PetscScalar *c,*c_q; 3795 3796 PetscFunctionBegin; 3797 m = A->rmap->n; 3798 n = A->cmap->n; 3799 p = B->cmap->n; 3800 a = sub_a->v; 3801 b = sub_b->a; 3802 c = sub_c->v; 3803 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3804 3805 ii = sub_b->i; 3806 idx = sub_b->j; 3807 for (i=0; i<n; i++) { 3808 q = ii[i+1] - ii[i]; 3809 while (q-->0) { 3810 c_q = c + m*(*idx); 3811 a_q = a + m*i; 3812 PetscKernelAXPY(c_q,*b,a_q,m); 3813 idx++; 3814 b++; 3815 } 3816 } 3817 PetscFunctionReturn(0); 3818 } 3819 3820 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3821 { 3822 PetscErrorCode ierr; 3823 PetscInt m=A->rmap->n,n=B->cmap->n; 3824 Mat Cmat; 3825 3826 PetscFunctionBegin; 3827 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); 3828 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3829 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3830 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3831 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3832 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3833 3834 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3835 3836 *C = Cmat; 3837 PetscFunctionReturn(0); 3838 } 3839 3840 /* ----------------------------------------------------------------*/ 3841 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3842 { 3843 PetscErrorCode ierr; 3844 3845 PetscFunctionBegin; 3846 if (scall == MAT_INITIAL_MATRIX) { 3847 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3848 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3849 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3850 } 3851 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3852 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3853 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3854 PetscFunctionReturn(0); 3855 } 3856 3857 3858 /*MC 3859 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3860 based on compressed sparse row format. 3861 3862 Options Database Keys: 3863 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3864 3865 Level: beginner 3866 3867 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3868 M*/ 3869 3870 /*MC 3871 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3872 3873 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3874 and MATMPIAIJ otherwise. As a result, for single process communicators, 3875 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3876 for communicators controlling multiple processes. It is recommended that you call both of 3877 the above preallocation routines for simplicity. 3878 3879 Options Database Keys: 3880 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3881 3882 Developer Notes: 3883 Subclasses include MATAIJCUSPARSE, MATAIJPERM, MATAIJMKL, MATAIJCRL, and also automatically switches over to use inodes when 3884 enough exist. 3885 3886 Level: beginner 3887 3888 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3889 M*/ 3890 3891 /*MC 3892 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3893 3894 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3895 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3896 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3897 for communicators controlling multiple processes. It is recommended that you call both of 3898 the above preallocation routines for simplicity. 3899 3900 Options Database Keys: 3901 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3902 3903 Level: beginner 3904 3905 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3906 M*/ 3907 3908 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3909 #if defined(PETSC_HAVE_ELEMENTAL) 3910 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3911 #endif 3912 #if defined(PETSC_HAVE_HYPRE) 3913 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3914 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3915 #endif 3916 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3917 3918 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3919 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3920 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3921 #endif 3922 3923 PETSC_EXTERN PetscErrorCode MatConvert_SeqAIJ_SeqSELL(Mat,MatType,MatReuse,Mat*); 3924 PETSC_INTERN PetscErrorCode MatConvert_XAIJ_IS(Mat,MatType,MatReuse,Mat*); 3925 PETSC_INTERN PetscErrorCode MatPtAP_IS_XAIJ(Mat,Mat,MatReuse,PetscReal,Mat*); 3926 3927 /*@C 3928 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 3929 3930 Not Collective 3931 3932 Input Parameter: 3933 . mat - a MATSEQAIJ matrix 3934 3935 Output Parameter: 3936 . array - pointer to the data 3937 3938 Level: intermediate 3939 3940 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3941 @*/ 3942 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3943 { 3944 PetscErrorCode ierr; 3945 3946 PetscFunctionBegin; 3947 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3948 PetscFunctionReturn(0); 3949 } 3950 3951 /*@C 3952 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3953 3954 Not Collective 3955 3956 Input Parameter: 3957 . mat - a MATSEQAIJ matrix 3958 3959 Output Parameter: 3960 . nz - the maximum number of nonzeros in any row 3961 3962 Level: intermediate 3963 3964 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3965 @*/ 3966 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3967 { 3968 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3969 3970 PetscFunctionBegin; 3971 *nz = aij->rmax; 3972 PetscFunctionReturn(0); 3973 } 3974 3975 /*@C 3976 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3977 3978 Not Collective 3979 3980 Input Parameters: 3981 . mat - a MATSEQAIJ matrix 3982 . array - pointer to the data 3983 3984 Level: intermediate 3985 3986 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3987 @*/ 3988 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3989 { 3990 PetscErrorCode ierr; 3991 3992 PetscFunctionBegin; 3993 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3994 PetscFunctionReturn(0); 3995 } 3996 3997 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3998 { 3999 Mat_SeqAIJ *b; 4000 PetscErrorCode ierr; 4001 PetscMPIInt size; 4002 4003 PetscFunctionBegin; 4004 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 4005 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 4006 4007 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 4008 4009 B->data = (void*)b; 4010 4011 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 4012 4013 b->row = 0; 4014 b->col = 0; 4015 b->icol = 0; 4016 b->reallocs = 0; 4017 b->ignorezeroentries = PETSC_FALSE; 4018 b->roworiented = PETSC_TRUE; 4019 b->nonew = 0; 4020 b->diag = 0; 4021 b->solve_work = 0; 4022 B->spptr = 0; 4023 b->saved_values = 0; 4024 b->idiag = 0; 4025 b->mdiag = 0; 4026 b->ssor_work = 0; 4027 b->omega = 1.0; 4028 b->fshift = 0.0; 4029 b->idiagvalid = PETSC_FALSE; 4030 b->ibdiagvalid = PETSC_FALSE; 4031 b->keepnonzeropattern = PETSC_FALSE; 4032 4033 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4034 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 4035 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 4036 4037 #if defined(PETSC_HAVE_MATLAB_ENGINE) 4038 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 4039 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 4040 #endif 4041 4042 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 4043 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 4044 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 4045 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 4046 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 4047 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4048 #if defined(PETSC_HAVE_MKL_SPARSE) 4049 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijmkl_C",MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4050 #endif 4051 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4052 #if defined(PETSC_HAVE_ELEMENTAL) 4053 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 4054 #endif 4055 #if defined(PETSC_HAVE_HYPRE) 4056 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 4057 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 4058 #endif 4059 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 4060 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsell_C",MatConvert_SeqAIJ_SeqSELL);CHKERRQ(ierr); 4061 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_is_C",MatConvert_XAIJ_IS);CHKERRQ(ierr); 4062 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4063 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 4064 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 4065 ierr = PetscObjectComposeFunction((PetscObject)B,"MatResetPreallocation_C",MatResetPreallocation_SeqAIJ);CHKERRQ(ierr); 4066 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 4067 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 4068 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 4069 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 4070 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 4071 ierr = PetscObjectComposeFunction((PetscObject)B,"MatPtAP_is_seqaij_C",MatPtAP_IS_XAIJ);CHKERRQ(ierr); 4072 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 4073 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 4074 ierr = MatSeqAIJSetTypeFromOptions(B);CHKERRQ(ierr); /* this allows changing the matrix subtype to say MATSEQAIJPERM */ 4075 PetscFunctionReturn(0); 4076 } 4077 4078 /* 4079 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4080 */ 4081 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4082 { 4083 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4084 PetscErrorCode ierr; 4085 PetscInt i,m = A->rmap->n; 4086 4087 PetscFunctionBegin; 4088 c = (Mat_SeqAIJ*)C->data; 4089 4090 C->factortype = A->factortype; 4091 c->row = 0; 4092 c->col = 0; 4093 c->icol = 0; 4094 c->reallocs = 0; 4095 4096 C->assembled = PETSC_TRUE; 4097 4098 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4099 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4100 4101 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4102 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4103 for (i=0; i<m; i++) { 4104 c->imax[i] = a->imax[i]; 4105 c->ilen[i] = a->ilen[i]; 4106 } 4107 4108 /* allocate the matrix space */ 4109 if (mallocmatspace) { 4110 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4111 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4112 4113 c->singlemalloc = PETSC_TRUE; 4114 4115 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4116 if (m > 0) { 4117 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4118 if (cpvalues == MAT_COPY_VALUES) { 4119 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4120 } else { 4121 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4122 } 4123 } 4124 } 4125 4126 c->ignorezeroentries = a->ignorezeroentries; 4127 c->roworiented = a->roworiented; 4128 c->nonew = a->nonew; 4129 if (a->diag) { 4130 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4131 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4132 for (i=0; i<m; i++) { 4133 c->diag[i] = a->diag[i]; 4134 } 4135 } else c->diag = 0; 4136 4137 c->solve_work = 0; 4138 c->saved_values = 0; 4139 c->idiag = 0; 4140 c->ssor_work = 0; 4141 c->keepnonzeropattern = a->keepnonzeropattern; 4142 c->free_a = PETSC_TRUE; 4143 c->free_ij = PETSC_TRUE; 4144 4145 c->rmax = a->rmax; 4146 c->nz = a->nz; 4147 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4148 C->preallocated = PETSC_TRUE; 4149 4150 c->compressedrow.use = a->compressedrow.use; 4151 c->compressedrow.nrows = a->compressedrow.nrows; 4152 if (a->compressedrow.use) { 4153 i = a->compressedrow.nrows; 4154 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4155 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4156 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4157 } else { 4158 c->compressedrow.use = PETSC_FALSE; 4159 c->compressedrow.i = NULL; 4160 c->compressedrow.rindex = NULL; 4161 } 4162 c->nonzerorowcnt = a->nonzerorowcnt; 4163 C->nonzerostate = A->nonzerostate; 4164 4165 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4166 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4167 PetscFunctionReturn(0); 4168 } 4169 4170 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4171 { 4172 PetscErrorCode ierr; 4173 4174 PetscFunctionBegin; 4175 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4176 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4177 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4178 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4179 } 4180 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4181 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4182 PetscFunctionReturn(0); 4183 } 4184 4185 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4186 { 4187 Mat_SeqAIJ *a; 4188 PetscErrorCode ierr; 4189 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4190 int fd; 4191 PetscMPIInt size; 4192 MPI_Comm comm; 4193 PetscInt bs = newMat->rmap->bs; 4194 4195 PetscFunctionBegin; 4196 /* force binary viewer to load .info file if it has not yet done so */ 4197 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4198 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4199 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4200 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4201 4202 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4203 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4204 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4205 if (bs < 0) bs = 1; 4206 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4207 4208 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4209 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4210 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4211 M = header[1]; N = header[2]; nz = header[3]; 4212 4213 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4214 4215 /* read in row lengths */ 4216 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4217 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4218 4219 /* check if sum of rowlengths is same as nz */ 4220 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4221 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); 4222 4223 /* set global size if not set already*/ 4224 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4225 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4226 } else { 4227 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4228 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4229 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4230 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4231 } 4232 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); 4233 } 4234 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4235 a = (Mat_SeqAIJ*)newMat->data; 4236 4237 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4238 4239 /* read in nonzero values */ 4240 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4241 4242 /* set matrix "i" values */ 4243 a->i[0] = 0; 4244 for (i=1; i<= M; i++) { 4245 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4246 a->ilen[i-1] = rowlengths[i-1]; 4247 } 4248 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4249 4250 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4251 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4252 PetscFunctionReturn(0); 4253 } 4254 4255 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4256 { 4257 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4258 PetscErrorCode ierr; 4259 #if defined(PETSC_USE_COMPLEX) 4260 PetscInt k; 4261 #endif 4262 4263 PetscFunctionBegin; 4264 /* If the matrix dimensions are not equal,or no of nonzeros */ 4265 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4266 *flg = PETSC_FALSE; 4267 PetscFunctionReturn(0); 4268 } 4269 4270 /* if the a->i are the same */ 4271 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4272 if (!*flg) PetscFunctionReturn(0); 4273 4274 /* if a->j are the same */ 4275 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4276 if (!*flg) PetscFunctionReturn(0); 4277 4278 /* if a->a are the same */ 4279 #if defined(PETSC_USE_COMPLEX) 4280 for (k=0; k<a->nz; k++) { 4281 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4282 *flg = PETSC_FALSE; 4283 PetscFunctionReturn(0); 4284 } 4285 } 4286 #else 4287 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4288 #endif 4289 PetscFunctionReturn(0); 4290 } 4291 4292 /*@ 4293 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4294 provided by the user. 4295 4296 Collective on MPI_Comm 4297 4298 Input Parameters: 4299 + comm - must be an MPI communicator of size 1 4300 . m - number of rows 4301 . n - number of columns 4302 . i - row indices 4303 . j - column indices 4304 - a - matrix values 4305 4306 Output Parameter: 4307 . mat - the matrix 4308 4309 Level: intermediate 4310 4311 Notes: 4312 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4313 once the matrix is destroyed and not before 4314 4315 You cannot set new nonzero locations into this matrix, that will generate an error. 4316 4317 The i and j indices are 0 based 4318 4319 The format which is used for the sparse matrix input, is equivalent to a 4320 row-major ordering.. i.e for the following matrix, the input data expected is 4321 as shown 4322 4323 $ 1 0 0 4324 $ 2 0 3 4325 $ 4 5 6 4326 $ 4327 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4328 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4329 $ v = {1,2,3,4,5,6} [size = 6] 4330 4331 4332 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4333 4334 @*/ 4335 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4336 { 4337 PetscErrorCode ierr; 4338 PetscInt ii; 4339 Mat_SeqAIJ *aij; 4340 #if defined(PETSC_USE_DEBUG) 4341 PetscInt jj; 4342 #endif 4343 4344 PetscFunctionBegin; 4345 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4346 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4347 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4348 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4349 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4350 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4351 aij = (Mat_SeqAIJ*)(*mat)->data; 4352 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4353 4354 aij->i = i; 4355 aij->j = j; 4356 aij->a = a; 4357 aij->singlemalloc = PETSC_FALSE; 4358 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4359 aij->free_a = PETSC_FALSE; 4360 aij->free_ij = PETSC_FALSE; 4361 4362 for (ii=0; ii<m; ii++) { 4363 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4364 #if defined(PETSC_USE_DEBUG) 4365 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]); 4366 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4367 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); 4368 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); 4369 } 4370 #endif 4371 } 4372 #if defined(PETSC_USE_DEBUG) 4373 for (ii=0; ii<aij->i[m]; ii++) { 4374 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4375 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]); 4376 } 4377 #endif 4378 4379 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4380 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4381 PetscFunctionReturn(0); 4382 } 4383 /*@C 4384 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4385 provided by the user. 4386 4387 Collective on MPI_Comm 4388 4389 Input Parameters: 4390 + comm - must be an MPI communicator of size 1 4391 . m - number of rows 4392 . n - number of columns 4393 . i - row indices 4394 . j - column indices 4395 . a - matrix values 4396 . nz - number of nonzeros 4397 - idx - 0 or 1 based 4398 4399 Output Parameter: 4400 . mat - the matrix 4401 4402 Level: intermediate 4403 4404 Notes: 4405 The i and j indices are 0 based 4406 4407 The format which is used for the sparse matrix input, is equivalent to a 4408 row-major ordering.. i.e for the following matrix, the input data expected is 4409 as shown: 4410 4411 1 0 0 4412 2 0 3 4413 4 5 6 4414 4415 i = {0,1,1,2,2,2} 4416 j = {0,0,2,0,1,2} 4417 v = {1,2,3,4,5,6} 4418 4419 4420 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4421 4422 @*/ 4423 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4424 { 4425 PetscErrorCode ierr; 4426 PetscInt ii, *nnz, one = 1,row,col; 4427 4428 4429 PetscFunctionBegin; 4430 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4431 for (ii = 0; ii < nz; ii++) { 4432 nnz[i[ii] - !!idx] += 1; 4433 } 4434 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4435 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4436 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4437 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4438 for (ii = 0; ii < nz; ii++) { 4439 if (idx) { 4440 row = i[ii] - 1; 4441 col = j[ii] - 1; 4442 } else { 4443 row = i[ii]; 4444 col = j[ii]; 4445 } 4446 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4447 } 4448 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4449 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4450 ierr = PetscFree(nnz);CHKERRQ(ierr); 4451 PetscFunctionReturn(0); 4452 } 4453 4454 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4455 { 4456 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4457 PetscErrorCode ierr; 4458 4459 PetscFunctionBegin; 4460 a->idiagvalid = PETSC_FALSE; 4461 a->ibdiagvalid = PETSC_FALSE; 4462 4463 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4464 PetscFunctionReturn(0); 4465 } 4466 4467 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4468 { 4469 PetscErrorCode ierr; 4470 PetscMPIInt size; 4471 4472 PetscFunctionBegin; 4473 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4474 if (size == 1) { 4475 if (scall == MAT_INITIAL_MATRIX) { 4476 ierr = MatDuplicate(inmat,MAT_COPY_VALUES,outmat);CHKERRQ(ierr); 4477 } else { 4478 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4479 } 4480 } else { 4481 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4482 } 4483 PetscFunctionReturn(0); 4484 } 4485 4486 /* 4487 Permute A into C's *local* index space using rowemb,colemb. 4488 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4489 of [0,m), colemb is in [0,n). 4490 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4491 */ 4492 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4493 { 4494 /* If making this function public, change the error returned in this function away from _PLIB. */ 4495 PetscErrorCode ierr; 4496 Mat_SeqAIJ *Baij; 4497 PetscBool seqaij; 4498 PetscInt m,n,*nz,i,j,count; 4499 PetscScalar v; 4500 const PetscInt *rowindices,*colindices; 4501 4502 PetscFunctionBegin; 4503 if (!B) PetscFunctionReturn(0); 4504 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4505 ierr = PetscObjectBaseTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4506 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4507 if (rowemb) { 4508 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4509 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); 4510 } else { 4511 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4512 } 4513 if (colemb) { 4514 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4515 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); 4516 } else { 4517 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4518 } 4519 4520 Baij = (Mat_SeqAIJ*)(B->data); 4521 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4522 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4523 for (i=0; i<B->rmap->n; i++) { 4524 nz[i] = Baij->i[i+1] - Baij->i[i]; 4525 } 4526 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4527 ierr = PetscFree(nz);CHKERRQ(ierr); 4528 } 4529 if (pattern == SUBSET_NONZERO_PATTERN) { 4530 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4531 } 4532 count = 0; 4533 rowindices = NULL; 4534 colindices = NULL; 4535 if (rowemb) { 4536 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4537 } 4538 if (colemb) { 4539 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4540 } 4541 for (i=0; i<B->rmap->n; i++) { 4542 PetscInt row; 4543 row = i; 4544 if (rowindices) row = rowindices[i]; 4545 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4546 PetscInt col; 4547 col = Baij->j[count]; 4548 if (colindices) col = colindices[col]; 4549 v = Baij->a[count]; 4550 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4551 ++count; 4552 } 4553 } 4554 /* FIXME: set C's nonzerostate correctly. */ 4555 /* Assembly for C is necessary. */ 4556 C->preallocated = PETSC_TRUE; 4557 C->assembled = PETSC_TRUE; 4558 C->was_assembled = PETSC_FALSE; 4559 PetscFunctionReturn(0); 4560 } 4561 4562 PetscFunctionList MatSeqAIJList = NULL; 4563 4564 /*@C 4565 MatSeqAIJSetType - Converts a MATSEQAIJ matrix to a subtype 4566 4567 Collective on Mat 4568 4569 Input Parameters: 4570 + mat - the matrix object 4571 - matype - matrix type 4572 4573 Options Database Key: 4574 . -mat_seqai_type <method> - for example seqaijcrl 4575 4576 4577 Level: intermediate 4578 4579 .keywords: Mat, MatType, set, method 4580 4581 .seealso: PCSetType(), VecSetType(), MatCreate(), MatType, Mat 4582 @*/ 4583 PetscErrorCode MatSeqAIJSetType(Mat mat, MatType matype) 4584 { 4585 PetscErrorCode ierr,(*r)(Mat,MatType,MatReuse,Mat*); 4586 PetscBool sametype; 4587 4588 PetscFunctionBegin; 4589 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 4590 ierr = PetscObjectTypeCompare((PetscObject)mat,matype,&sametype);CHKERRQ(ierr); 4591 if (sametype) PetscFunctionReturn(0); 4592 4593 ierr = PetscFunctionListFind(MatSeqAIJList,matype,&r);CHKERRQ(ierr); 4594 if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown Mat type given: %s",matype); 4595 ierr = (*r)(mat,matype,MAT_INPLACE_MATRIX,&mat);CHKERRQ(ierr); 4596 PetscFunctionReturn(0); 4597 } 4598 4599 4600 /*@C 4601 MatSeqAIJRegister - - Adds a new sub-matrix type for sequential AIJ matrices 4602 4603 Not Collective 4604 4605 Input Parameters: 4606 + name - name of a new user-defined matrix type, for example MATSEQAIJCRL 4607 - function - routine to convert to subtype 4608 4609 Notes: 4610 MatSeqAIJRegister() may be called multiple times to add several user-defined solvers. 4611 4612 4613 Then, your matrix can be chosen with the procedural interface at runtime via the option 4614 $ -mat_seqaij_type my_mat 4615 4616 Level: advanced 4617 4618 .keywords: Mat, register 4619 4620 .seealso: MatSeqAIJRegisterAll() 4621 4622 4623 Level: advanced 4624 @*/ 4625 PetscErrorCode MatSeqAIJRegister(const char sname[],PetscErrorCode (*function)(Mat,MatType,MatReuse,Mat *)) 4626 { 4627 PetscErrorCode ierr; 4628 4629 PetscFunctionBegin; 4630 ierr = PetscFunctionListAdd(&MatSeqAIJList,sname,function);CHKERRQ(ierr); 4631 PetscFunctionReturn(0); 4632 } 4633 4634 PetscBool MatSeqAIJRegisterAllCalled = PETSC_FALSE; 4635 4636 /*@C 4637 MatSeqAIJRegisterAll - Registers all of the matrix subtypes of SeqAIJ 4638 4639 Not Collective 4640 4641 Level: advanced 4642 4643 Developers Note: CUSP and CUSPARSE do not yet support the MatConvert_SeqAIJ..() paradigm and thus cannot be registered here 4644 4645 .keywords: KSP, register, all 4646 4647 .seealso: MatRegisterAll(), MatSeqAIJRegister() 4648 @*/ 4649 PetscErrorCode MatSeqAIJRegisterAll(void) 4650 { 4651 PetscErrorCode ierr; 4652 4653 PetscFunctionBegin; 4654 if (MatSeqAIJRegisterAllCalled) PetscFunctionReturn(0); 4655 MatSeqAIJRegisterAllCalled = PETSC_TRUE; 4656 4657 ierr = MatSeqAIJRegister(MATSEQAIJCRL, MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 4658 ierr = MatSeqAIJRegister(MATSEQAIJPERM, MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 4659 #if defined(PETSC_HAVE_MKL_SPARSE) 4660 ierr = MatSeqAIJRegister(MATSEQAIJMKL, MatConvert_SeqAIJ_SeqAIJMKL);CHKERRQ(ierr); 4661 #endif 4662 #if defined(PETSC_HAVE_VIENNACL) && defined(PETSC_HAVE_VIENNACL_NO_CUDA) 4663 ierr = MatSeqAIJRegister(MATMPIAIJVIENNACL, MatConvert_SeqAIJ_SeqAIJViennaCL);CHKERRQ(ierr); 4664 #endif 4665 PetscFunctionReturn(0); 4666 } 4667 4668 /* 4669 Special version for direct calls from Fortran 4670 */ 4671 #include <petsc/private/fortranimpl.h> 4672 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4673 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4674 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4675 #define matsetvaluesseqaij_ matsetvaluesseqaij 4676 #endif 4677 4678 /* Change these macros so can be used in void function */ 4679 #undef CHKERRQ 4680 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4681 #undef SETERRQ2 4682 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4683 #undef SETERRQ3 4684 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4685 4686 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) 4687 { 4688 Mat A = *AA; 4689 PetscInt m = *mm, n = *nn; 4690 InsertMode is = *isis; 4691 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4692 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4693 PetscInt *imax,*ai,*ailen; 4694 PetscErrorCode ierr; 4695 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4696 MatScalar *ap,value,*aa; 4697 PetscBool ignorezeroentries = a->ignorezeroentries; 4698 PetscBool roworiented = a->roworiented; 4699 4700 PetscFunctionBegin; 4701 MatCheckPreallocated(A,1); 4702 imax = a->imax; 4703 ai = a->i; 4704 ailen = a->ilen; 4705 aj = a->j; 4706 aa = a->a; 4707 4708 for (k=0; k<m; k++) { /* loop over added rows */ 4709 row = im[k]; 4710 if (row < 0) continue; 4711 #if defined(PETSC_USE_DEBUG) 4712 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4713 #endif 4714 rp = aj + ai[row]; ap = aa + ai[row]; 4715 rmax = imax[row]; nrow = ailen[row]; 4716 low = 0; 4717 high = nrow; 4718 for (l=0; l<n; l++) { /* loop over added columns */ 4719 if (in[l] < 0) continue; 4720 #if defined(PETSC_USE_DEBUG) 4721 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4722 #endif 4723 col = in[l]; 4724 if (roworiented) value = v[l + k*n]; 4725 else value = v[k + l*m]; 4726 4727 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4728 4729 if (col <= lastcol) low = 0; 4730 else high = nrow; 4731 lastcol = col; 4732 while (high-low > 5) { 4733 t = (low+high)/2; 4734 if (rp[t] > col) high = t; 4735 else low = t; 4736 } 4737 for (i=low; i<high; i++) { 4738 if (rp[i] > col) break; 4739 if (rp[i] == col) { 4740 if (is == ADD_VALUES) ap[i] += value; 4741 else ap[i] = value; 4742 goto noinsert; 4743 } 4744 } 4745 if (value == 0.0 && ignorezeroentries) goto noinsert; 4746 if (nonew == 1) goto noinsert; 4747 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4748 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4749 N = nrow++ - 1; a->nz++; high++; 4750 /* shift up all the later entries in this row */ 4751 for (ii=N; ii>=i; ii--) { 4752 rp[ii+1] = rp[ii]; 4753 ap[ii+1] = ap[ii]; 4754 } 4755 rp[i] = col; 4756 ap[i] = value; 4757 A->nonzerostate++; 4758 noinsert:; 4759 low = i + 1; 4760 } 4761 ailen[row] = nrow; 4762 } 4763 PetscFunctionReturnVoid(); 4764 } 4765