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