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