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