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 ierr = PetscObjectStateIncrease((PetscObject)B);CHKERRQ(ierr); 2643 } else { 2644 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2645 } 2646 PetscFunctionReturn(0); 2647 } 2648 2649 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2650 { 2651 PetscErrorCode ierr; 2652 2653 PetscFunctionBegin; 2654 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2655 PetscFunctionReturn(0); 2656 } 2657 2658 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2659 { 2660 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2661 2662 PetscFunctionBegin; 2663 *array = a->a; 2664 PetscFunctionReturn(0); 2665 } 2666 2667 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2668 { 2669 PetscFunctionBegin; 2670 PetscFunctionReturn(0); 2671 } 2672 2673 /* 2674 Computes the number of nonzeros per row needed for preallocation when X and Y 2675 have different nonzero structure. 2676 */ 2677 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2678 { 2679 PetscInt i,j,k,nzx,nzy; 2680 2681 PetscFunctionBegin; 2682 /* Set the number of nonzeros in the new matrix */ 2683 for (i=0; i<m; i++) { 2684 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2685 nzx = xi[i+1] - xi[i]; 2686 nzy = yi[i+1] - yi[i]; 2687 nnz[i] = 0; 2688 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2689 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2690 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2691 nnz[i]++; 2692 } 2693 for (; k<nzy; k++) nnz[i]++; 2694 } 2695 PetscFunctionReturn(0); 2696 } 2697 2698 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2699 { 2700 PetscInt m = Y->rmap->N; 2701 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2702 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2703 PetscErrorCode ierr; 2704 2705 PetscFunctionBegin; 2706 /* Set the number of nonzeros in the new matrix */ 2707 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2708 PetscFunctionReturn(0); 2709 } 2710 2711 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2712 { 2713 PetscErrorCode ierr; 2714 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2715 PetscBLASInt one=1,bnz; 2716 2717 PetscFunctionBegin; 2718 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2719 if (str == SAME_NONZERO_PATTERN) { 2720 PetscScalar alpha = a; 2721 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2722 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2723 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2724 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2725 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2726 } else { 2727 Mat B; 2728 PetscInt *nnz; 2729 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2730 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2731 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2732 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2733 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2734 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2735 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2736 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2737 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2738 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2739 ierr = PetscFree(nnz);CHKERRQ(ierr); 2740 } 2741 PetscFunctionReturn(0); 2742 } 2743 2744 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2745 { 2746 #if defined(PETSC_USE_COMPLEX) 2747 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2748 PetscInt i,nz; 2749 PetscScalar *a; 2750 2751 PetscFunctionBegin; 2752 nz = aij->nz; 2753 a = aij->a; 2754 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2755 #else 2756 PetscFunctionBegin; 2757 #endif 2758 PetscFunctionReturn(0); 2759 } 2760 2761 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2762 { 2763 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2764 PetscErrorCode ierr; 2765 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2766 PetscReal atmp; 2767 PetscScalar *x; 2768 MatScalar *aa; 2769 2770 PetscFunctionBegin; 2771 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2772 aa = a->a; 2773 ai = a->i; 2774 aj = a->j; 2775 2776 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2777 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2778 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2779 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2780 for (i=0; i<m; i++) { 2781 ncols = ai[1] - ai[0]; ai++; 2782 x[i] = 0.0; 2783 for (j=0; j<ncols; j++) { 2784 atmp = PetscAbsScalar(*aa); 2785 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2786 aa++; aj++; 2787 } 2788 } 2789 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2790 PetscFunctionReturn(0); 2791 } 2792 2793 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2794 { 2795 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2796 PetscErrorCode ierr; 2797 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2798 PetscScalar *x; 2799 MatScalar *aa; 2800 2801 PetscFunctionBegin; 2802 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2803 aa = a->a; 2804 ai = a->i; 2805 aj = a->j; 2806 2807 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2808 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2809 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2810 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2811 for (i=0; i<m; i++) { 2812 ncols = ai[1] - ai[0]; ai++; 2813 if (ncols == A->cmap->n) { /* row is dense */ 2814 x[i] = *aa; if (idx) idx[i] = 0; 2815 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2816 x[i] = 0.0; 2817 if (idx) { 2818 idx[i] = 0; /* in case ncols is zero */ 2819 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2820 if (aj[j] > j) { 2821 idx[i] = j; 2822 break; 2823 } 2824 } 2825 } 2826 } 2827 for (j=0; j<ncols; j++) { 2828 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2829 aa++; aj++; 2830 } 2831 } 2832 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2833 PetscFunctionReturn(0); 2834 } 2835 2836 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2837 { 2838 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2839 PetscErrorCode ierr; 2840 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2841 PetscReal atmp; 2842 PetscScalar *x; 2843 MatScalar *aa; 2844 2845 PetscFunctionBegin; 2846 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2847 aa = a->a; 2848 ai = a->i; 2849 aj = a->j; 2850 2851 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2852 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2853 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2854 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); 2855 for (i=0; i<m; i++) { 2856 ncols = ai[1] - ai[0]; ai++; 2857 if (ncols) { 2858 /* Get first nonzero */ 2859 for (j = 0; j < ncols; j++) { 2860 atmp = PetscAbsScalar(aa[j]); 2861 if (atmp > 1.0e-12) { 2862 x[i] = atmp; 2863 if (idx) idx[i] = aj[j]; 2864 break; 2865 } 2866 } 2867 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2868 } else { 2869 x[i] = 0.0; if (idx) idx[i] = 0; 2870 } 2871 for (j = 0; j < ncols; j++) { 2872 atmp = PetscAbsScalar(*aa); 2873 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2874 aa++; aj++; 2875 } 2876 } 2877 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2878 PetscFunctionReturn(0); 2879 } 2880 2881 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2882 { 2883 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2884 PetscErrorCode ierr; 2885 PetscInt i,j,m = A->rmap->n,ncols,n; 2886 const PetscInt *ai,*aj; 2887 PetscScalar *x; 2888 const MatScalar *aa; 2889 2890 PetscFunctionBegin; 2891 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2892 aa = a->a; 2893 ai = a->i; 2894 aj = a->j; 2895 2896 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2897 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2898 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2899 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2900 for (i=0; i<m; i++) { 2901 ncols = ai[1] - ai[0]; ai++; 2902 if (ncols == A->cmap->n) { /* row is dense */ 2903 x[i] = *aa; if (idx) idx[i] = 0; 2904 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2905 x[i] = 0.0; 2906 if (idx) { /* find first implicit 0.0 in the row */ 2907 idx[i] = 0; /* in case ncols is zero */ 2908 for (j=0; j<ncols; j++) { 2909 if (aj[j] > j) { 2910 idx[i] = j; 2911 break; 2912 } 2913 } 2914 } 2915 } 2916 for (j=0; j<ncols; j++) { 2917 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2918 aa++; aj++; 2919 } 2920 } 2921 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2922 PetscFunctionReturn(0); 2923 } 2924 2925 #include <petscblaslapack.h> 2926 #include <petsc/private/kernels/blockinvert.h> 2927 2928 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2929 { 2930 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2931 PetscErrorCode ierr; 2932 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2933 MatScalar *diag,work[25],*v_work; 2934 PetscReal shift = 0.0; 2935 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2936 2937 PetscFunctionBegin; 2938 allowzeropivot = PetscNot(A->erroriffailure); 2939 if (a->ibdiagvalid) { 2940 if (values) *values = a->ibdiag; 2941 PetscFunctionReturn(0); 2942 } 2943 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2944 if (!a->ibdiag) { 2945 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2946 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2947 } 2948 diag = a->ibdiag; 2949 if (values) *values = a->ibdiag; 2950 /* factor and invert each block */ 2951 switch (bs) { 2952 case 1: 2953 for (i=0; i<mbs; i++) { 2954 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2955 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 2956 if (allowzeropivot) { 2957 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2958 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 2959 A->factorerror_zeropivot_row = i; 2960 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 2961 } 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); 2962 } 2963 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2964 } 2965 break; 2966 case 2: 2967 for (i=0; i<mbs; i++) { 2968 ij[0] = 2*i; ij[1] = 2*i + 1; 2969 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2970 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2971 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2972 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2973 diag += 4; 2974 } 2975 break; 2976 case 3: 2977 for (i=0; i<mbs; i++) { 2978 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2979 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2980 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2981 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2982 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2983 diag += 9; 2984 } 2985 break; 2986 case 4: 2987 for (i=0; i<mbs; i++) { 2988 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2989 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2990 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2991 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2992 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2993 diag += 16; 2994 } 2995 break; 2996 case 5: 2997 for (i=0; i<mbs; i++) { 2998 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2999 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 3000 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3001 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3002 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 3003 diag += 25; 3004 } 3005 break; 3006 case 6: 3007 for (i=0; i<mbs; i++) { 3008 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; 3009 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 3010 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3011 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3012 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 3013 diag += 36; 3014 } 3015 break; 3016 case 7: 3017 for (i=0; i<mbs; i++) { 3018 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; 3019 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 3020 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3021 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3022 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 3023 diag += 49; 3024 } 3025 break; 3026 default: 3027 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 3028 for (i=0; i<mbs; i++) { 3029 for (j=0; j<bs; j++) { 3030 IJ[j] = bs*i + j; 3031 } 3032 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3033 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 3034 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 3035 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3036 diag += bs2; 3037 } 3038 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3039 } 3040 a->ibdiagvalid = PETSC_TRUE; 3041 PetscFunctionReturn(0); 3042 } 3043 3044 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3045 { 3046 PetscErrorCode ierr; 3047 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3048 PetscScalar a; 3049 PetscInt m,n,i,j,col; 3050 3051 PetscFunctionBegin; 3052 if (!x->assembled) { 3053 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3054 for (i=0; i<m; i++) { 3055 for (j=0; j<aij->imax[i]; j++) { 3056 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3057 col = (PetscInt)(n*PetscRealPart(a)); 3058 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3059 } 3060 } 3061 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3062 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3063 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3064 PetscFunctionReturn(0); 3065 } 3066 3067 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3068 { 3069 PetscErrorCode ierr; 3070 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3071 3072 PetscFunctionBegin; 3073 if (!Y->preallocated || !aij->nz) { 3074 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3075 } 3076 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3077 PetscFunctionReturn(0); 3078 } 3079 3080 /* -------------------------------------------------------------------*/ 3081 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3082 MatGetRow_SeqAIJ, 3083 MatRestoreRow_SeqAIJ, 3084 MatMult_SeqAIJ, 3085 /* 4*/ MatMultAdd_SeqAIJ, 3086 MatMultTranspose_SeqAIJ, 3087 MatMultTransposeAdd_SeqAIJ, 3088 0, 3089 0, 3090 0, 3091 /* 10*/ 0, 3092 MatLUFactor_SeqAIJ, 3093 0, 3094 MatSOR_SeqAIJ, 3095 MatTranspose_SeqAIJ, 3096 /*1 5*/ MatGetInfo_SeqAIJ, 3097 MatEqual_SeqAIJ, 3098 MatGetDiagonal_SeqAIJ, 3099 MatDiagonalScale_SeqAIJ, 3100 MatNorm_SeqAIJ, 3101 /* 20*/ 0, 3102 MatAssemblyEnd_SeqAIJ, 3103 MatSetOption_SeqAIJ, 3104 MatZeroEntries_SeqAIJ, 3105 /* 24*/ MatZeroRows_SeqAIJ, 3106 0, 3107 0, 3108 0, 3109 0, 3110 /* 29*/ MatSetUp_SeqAIJ, 3111 0, 3112 0, 3113 0, 3114 0, 3115 /* 34*/ MatDuplicate_SeqAIJ, 3116 0, 3117 0, 3118 MatILUFactor_SeqAIJ, 3119 0, 3120 /* 39*/ MatAXPY_SeqAIJ, 3121 MatCreateSubMatrices_SeqAIJ, 3122 MatIncreaseOverlap_SeqAIJ, 3123 MatGetValues_SeqAIJ, 3124 MatCopy_SeqAIJ, 3125 /* 44*/ MatGetRowMax_SeqAIJ, 3126 MatScale_SeqAIJ, 3127 MatShift_SeqAIJ, 3128 MatDiagonalSet_SeqAIJ, 3129 MatZeroRowsColumns_SeqAIJ, 3130 /* 49*/ MatSetRandom_SeqAIJ, 3131 MatGetRowIJ_SeqAIJ, 3132 MatRestoreRowIJ_SeqAIJ, 3133 MatGetColumnIJ_SeqAIJ, 3134 MatRestoreColumnIJ_SeqAIJ, 3135 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3136 0, 3137 0, 3138 MatPermute_SeqAIJ, 3139 0, 3140 /* 59*/ 0, 3141 MatDestroy_SeqAIJ, 3142 MatView_SeqAIJ, 3143 0, 3144 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3145 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3146 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3147 0, 3148 0, 3149 0, 3150 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3151 MatGetRowMinAbs_SeqAIJ, 3152 0, 3153 0, 3154 0, 3155 /* 74*/ 0, 3156 MatFDColoringApply_AIJ, 3157 0, 3158 0, 3159 0, 3160 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3161 0, 3162 0, 3163 0, 3164 MatLoad_SeqAIJ, 3165 /* 84*/ MatIsSymmetric_SeqAIJ, 3166 MatIsHermitian_SeqAIJ, 3167 0, 3168 0, 3169 0, 3170 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3171 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3172 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3173 MatPtAP_SeqAIJ_SeqAIJ, 3174 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3175 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3176 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3177 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3178 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3179 0, 3180 /* 99*/ 0, 3181 0, 3182 0, 3183 MatConjugate_SeqAIJ, 3184 0, 3185 /*104*/ MatSetValuesRow_SeqAIJ, 3186 MatRealPart_SeqAIJ, 3187 MatImaginaryPart_SeqAIJ, 3188 0, 3189 0, 3190 /*109*/ MatMatSolve_SeqAIJ, 3191 0, 3192 MatGetRowMin_SeqAIJ, 3193 0, 3194 MatMissingDiagonal_SeqAIJ, 3195 /*114*/ 0, 3196 0, 3197 0, 3198 0, 3199 0, 3200 /*119*/ 0, 3201 0, 3202 0, 3203 0, 3204 MatGetMultiProcBlock_SeqAIJ, 3205 /*124*/ MatFindNonzeroRows_SeqAIJ, 3206 MatGetColumnNorms_SeqAIJ, 3207 MatInvertBlockDiagonal_SeqAIJ, 3208 0, 3209 0, 3210 /*129*/ 0, 3211 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3212 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3213 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3214 MatTransposeColoringCreate_SeqAIJ, 3215 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3216 MatTransColoringApplyDenToSp_SeqAIJ, 3217 MatRARt_SeqAIJ_SeqAIJ, 3218 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3219 MatRARtNumeric_SeqAIJ_SeqAIJ, 3220 /*139*/0, 3221 0, 3222 0, 3223 MatFDColoringSetUp_SeqXAIJ, 3224 MatFindOffBlockDiagonalEntries_SeqAIJ, 3225 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ 3226 }; 3227 3228 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3229 { 3230 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3231 PetscInt i,nz,n; 3232 3233 PetscFunctionBegin; 3234 nz = aij->maxnz; 3235 n = mat->rmap->n; 3236 for (i=0; i<nz; i++) { 3237 aij->j[i] = indices[i]; 3238 } 3239 aij->nz = nz; 3240 for (i=0; i<n; i++) { 3241 aij->ilen[i] = aij->imax[i]; 3242 } 3243 PetscFunctionReturn(0); 3244 } 3245 3246 /*@ 3247 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3248 in the matrix. 3249 3250 Input Parameters: 3251 + mat - the SeqAIJ matrix 3252 - indices - the column indices 3253 3254 Level: advanced 3255 3256 Notes: 3257 This can be called if you have precomputed the nonzero structure of the 3258 matrix and want to provide it to the matrix object to improve the performance 3259 of the MatSetValues() operation. 3260 3261 You MUST have set the correct numbers of nonzeros per row in the call to 3262 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3263 3264 MUST be called before any calls to MatSetValues(); 3265 3266 The indices should start with zero, not one. 3267 3268 @*/ 3269 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3270 { 3271 PetscErrorCode ierr; 3272 3273 PetscFunctionBegin; 3274 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3275 PetscValidPointer(indices,2); 3276 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3277 PetscFunctionReturn(0); 3278 } 3279 3280 /* ----------------------------------------------------------------------------------------*/ 3281 3282 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3283 { 3284 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3285 PetscErrorCode ierr; 3286 size_t nz = aij->i[mat->rmap->n]; 3287 3288 PetscFunctionBegin; 3289 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3290 3291 /* allocate space for values if not already there */ 3292 if (!aij->saved_values) { 3293 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3294 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3295 } 3296 3297 /* copy values over */ 3298 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3299 PetscFunctionReturn(0); 3300 } 3301 3302 /*@ 3303 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3304 example, reuse of the linear part of a Jacobian, while recomputing the 3305 nonlinear portion. 3306 3307 Collect on Mat 3308 3309 Input Parameters: 3310 . mat - the matrix (currently only AIJ matrices support this option) 3311 3312 Level: advanced 3313 3314 Common Usage, with SNESSolve(): 3315 $ Create Jacobian matrix 3316 $ Set linear terms into matrix 3317 $ Apply boundary conditions to matrix, at this time matrix must have 3318 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3319 $ boundary conditions again will not change the nonzero structure 3320 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3321 $ ierr = MatStoreValues(mat); 3322 $ Call SNESSetJacobian() with matrix 3323 $ In your Jacobian routine 3324 $ ierr = MatRetrieveValues(mat); 3325 $ Set nonlinear terms in matrix 3326 3327 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3328 $ // build linear portion of Jacobian 3329 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3330 $ ierr = MatStoreValues(mat); 3331 $ loop over nonlinear iterations 3332 $ ierr = MatRetrieveValues(mat); 3333 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3334 $ // call MatAssemblyBegin/End() on matrix 3335 $ Solve linear system with Jacobian 3336 $ endloop 3337 3338 Notes: 3339 Matrix must already be assemblied before calling this routine 3340 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3341 calling this routine. 3342 3343 When this is called multiple times it overwrites the previous set of stored values 3344 and does not allocated additional space. 3345 3346 .seealso: MatRetrieveValues() 3347 3348 @*/ 3349 PetscErrorCode MatStoreValues(Mat mat) 3350 { 3351 PetscErrorCode ierr; 3352 3353 PetscFunctionBegin; 3354 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3355 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3356 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3357 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3358 PetscFunctionReturn(0); 3359 } 3360 3361 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3362 { 3363 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3364 PetscErrorCode ierr; 3365 PetscInt nz = aij->i[mat->rmap->n]; 3366 3367 PetscFunctionBegin; 3368 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3369 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3370 /* copy values over */ 3371 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3372 PetscFunctionReturn(0); 3373 } 3374 3375 /*@ 3376 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3377 example, reuse of the linear part of a Jacobian, while recomputing the 3378 nonlinear portion. 3379 3380 Collect on Mat 3381 3382 Input Parameters: 3383 . mat - the matrix (currently only AIJ matrices support this option) 3384 3385 Level: advanced 3386 3387 .seealso: MatStoreValues() 3388 3389 @*/ 3390 PetscErrorCode MatRetrieveValues(Mat mat) 3391 { 3392 PetscErrorCode ierr; 3393 3394 PetscFunctionBegin; 3395 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3396 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3397 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3398 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3399 PetscFunctionReturn(0); 3400 } 3401 3402 3403 /* --------------------------------------------------------------------------------*/ 3404 /*@C 3405 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3406 (the default parallel PETSc format). For good matrix assembly performance 3407 the user should preallocate the matrix storage by setting the parameter nz 3408 (or the array nnz). By setting these parameters accurately, performance 3409 during matrix assembly can be increased by more than a factor of 50. 3410 3411 Collective on MPI_Comm 3412 3413 Input Parameters: 3414 + comm - MPI communicator, set to PETSC_COMM_SELF 3415 . m - number of rows 3416 . n - number of columns 3417 . nz - number of nonzeros per row (same for all rows) 3418 - nnz - array containing the number of nonzeros in the various rows 3419 (possibly different for each row) or NULL 3420 3421 Output Parameter: 3422 . A - the matrix 3423 3424 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3425 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3426 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3427 3428 Notes: 3429 If nnz is given then nz is ignored 3430 3431 The AIJ format (also called the Yale sparse matrix format or 3432 compressed row storage), is fully compatible with standard Fortran 77 3433 storage. That is, the stored row and column indices can begin at 3434 either one (as in Fortran) or zero. See the users' manual for details. 3435 3436 Specify the preallocated storage with either nz or nnz (not both). 3437 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3438 allocation. For large problems you MUST preallocate memory or you 3439 will get TERRIBLE performance, see the users' manual chapter on matrices. 3440 3441 By default, this format uses inodes (identical nodes) when possible, to 3442 improve numerical efficiency of matrix-vector products and solves. We 3443 search for consecutive rows with the same nonzero structure, thereby 3444 reusing matrix information to achieve increased efficiency. 3445 3446 Options Database Keys: 3447 + -mat_no_inode - Do not use inodes 3448 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3449 3450 Level: intermediate 3451 3452 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3453 3454 @*/ 3455 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3456 { 3457 PetscErrorCode ierr; 3458 3459 PetscFunctionBegin; 3460 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3461 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3462 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3463 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3464 PetscFunctionReturn(0); 3465 } 3466 3467 /*@C 3468 MatSeqAIJSetPreallocation - For good matrix assembly performance 3469 the user should preallocate the matrix storage by setting the parameter nz 3470 (or the array nnz). By setting these parameters accurately, performance 3471 during matrix assembly can be increased by more than a factor of 50. 3472 3473 Collective on MPI_Comm 3474 3475 Input Parameters: 3476 + B - The matrix 3477 . nz - number of nonzeros per row (same for all rows) 3478 - nnz - array containing the number of nonzeros in the various rows 3479 (possibly different for each row) or NULL 3480 3481 Notes: 3482 If nnz is given then nz is ignored 3483 3484 The AIJ format (also called the Yale sparse matrix format or 3485 compressed row storage), is fully compatible with standard Fortran 77 3486 storage. That is, the stored row and column indices can begin at 3487 either one (as in Fortran) or zero. See the users' manual for details. 3488 3489 Specify the preallocated storage with either nz or nnz (not both). 3490 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3491 allocation. For large problems you MUST preallocate memory or you 3492 will get TERRIBLE performance, see the users' manual chapter on matrices. 3493 3494 You can call MatGetInfo() to get information on how effective the preallocation was; 3495 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3496 You can also run with the option -info and look for messages with the string 3497 malloc in them to see if additional memory allocation was needed. 3498 3499 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3500 entries or columns indices 3501 3502 By default, this format uses inodes (identical nodes) when possible, to 3503 improve numerical efficiency of matrix-vector products and solves. We 3504 search for consecutive rows with the same nonzero structure, thereby 3505 reusing matrix information to achieve increased efficiency. 3506 3507 Options Database Keys: 3508 + -mat_no_inode - Do not use inodes 3509 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3510 - -mat_aij_oneindex - Internally use indexing starting at 1 3511 rather than 0. Note that when calling MatSetValues(), 3512 the user still MUST index entries starting at 0! 3513 3514 Level: intermediate 3515 3516 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3517 3518 @*/ 3519 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3520 { 3521 PetscErrorCode ierr; 3522 3523 PetscFunctionBegin; 3524 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3525 PetscValidType(B,1); 3526 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3527 PetscFunctionReturn(0); 3528 } 3529 3530 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3531 { 3532 Mat_SeqAIJ *b; 3533 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3534 PetscErrorCode ierr; 3535 PetscInt i; 3536 3537 PetscFunctionBegin; 3538 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3539 if (nz == MAT_SKIP_ALLOCATION) { 3540 skipallocation = PETSC_TRUE; 3541 nz = 0; 3542 } 3543 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3544 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3545 3546 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3547 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3548 if (nnz) { 3549 for (i=0; i<B->rmap->n; i++) { 3550 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]); 3551 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); 3552 } 3553 } 3554 3555 B->preallocated = PETSC_TRUE; 3556 3557 b = (Mat_SeqAIJ*)B->data; 3558 3559 if (!skipallocation) { 3560 if (!b->imax) { 3561 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3562 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3563 } 3564 if (!nnz) { 3565 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3566 else if (nz < 0) nz = 1; 3567 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3568 nz = nz*B->rmap->n; 3569 } else { 3570 nz = 0; 3571 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3572 } 3573 /* b->ilen will count nonzeros in each row so far. */ 3574 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3575 3576 /* allocate the matrix space */ 3577 /* FIXME: should B's old memory be unlogged? */ 3578 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3579 if (B->structure_only) { 3580 ierr = PetscMalloc1(nz,&b->j);CHKERRQ(ierr); 3581 ierr = PetscMalloc1(B->rmap->n+1,&b->i);CHKERRQ(ierr); 3582 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3583 } else { 3584 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3585 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3586 } 3587 b->i[0] = 0; 3588 for (i=1; i<B->rmap->n+1; i++) { 3589 b->i[i] = b->i[i-1] + b->imax[i-1]; 3590 } 3591 if (B->structure_only) { 3592 b->singlemalloc = PETSC_FALSE; 3593 b->free_a = PETSC_FALSE; 3594 } else { 3595 b->singlemalloc = PETSC_TRUE; 3596 b->free_a = PETSC_TRUE; 3597 } 3598 b->free_ij = PETSC_TRUE; 3599 } else { 3600 b->free_a = PETSC_FALSE; 3601 b->free_ij = PETSC_FALSE; 3602 } 3603 3604 b->nz = 0; 3605 b->maxnz = nz; 3606 B->info.nz_unneeded = (double)b->maxnz; 3607 if (realalloc) { 3608 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3609 } 3610 B->was_assembled = PETSC_FALSE; 3611 B->assembled = PETSC_FALSE; 3612 PetscFunctionReturn(0); 3613 } 3614 3615 /*@ 3616 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3617 3618 Input Parameters: 3619 + B - the matrix 3620 . i - the indices into j for the start of each row (starts with zero) 3621 . j - the column indices for each row (starts with zero) these must be sorted for each row 3622 - v - optional values in the matrix 3623 3624 Level: developer 3625 3626 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3627 3628 .keywords: matrix, aij, compressed row, sparse, sequential 3629 3630 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3631 @*/ 3632 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3633 { 3634 PetscErrorCode ierr; 3635 3636 PetscFunctionBegin; 3637 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3638 PetscValidType(B,1); 3639 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3640 PetscFunctionReturn(0); 3641 } 3642 3643 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3644 { 3645 PetscInt i; 3646 PetscInt m,n; 3647 PetscInt nz; 3648 PetscInt *nnz, nz_max = 0; 3649 PetscScalar *values; 3650 PetscErrorCode ierr; 3651 3652 PetscFunctionBegin; 3653 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3654 3655 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3656 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3657 3658 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3659 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3660 for (i = 0; i < m; i++) { 3661 nz = Ii[i+1]- Ii[i]; 3662 nz_max = PetscMax(nz_max, nz); 3663 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3664 nnz[i] = nz; 3665 } 3666 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3667 ierr = PetscFree(nnz);CHKERRQ(ierr); 3668 3669 if (v) { 3670 values = (PetscScalar*) v; 3671 } else { 3672 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3673 } 3674 3675 for (i = 0; i < m; i++) { 3676 nz = Ii[i+1] - Ii[i]; 3677 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3678 } 3679 3680 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3681 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3682 3683 if (!v) { 3684 ierr = PetscFree(values);CHKERRQ(ierr); 3685 } 3686 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3687 PetscFunctionReturn(0); 3688 } 3689 3690 #include <../src/mat/impls/dense/seq/dense.h> 3691 #include <petsc/private/kernels/petscaxpy.h> 3692 3693 /* 3694 Computes (B'*A')' since computing B*A directly is untenable 3695 3696 n p p 3697 ( ) ( ) ( ) 3698 m ( A ) * n ( B ) = m ( C ) 3699 ( ) ( ) ( ) 3700 3701 */ 3702 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3703 { 3704 PetscErrorCode ierr; 3705 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3706 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3707 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3708 PetscInt i,n,m,q,p; 3709 const PetscInt *ii,*idx; 3710 const PetscScalar *b,*a,*a_q; 3711 PetscScalar *c,*c_q; 3712 3713 PetscFunctionBegin; 3714 m = A->rmap->n; 3715 n = A->cmap->n; 3716 p = B->cmap->n; 3717 a = sub_a->v; 3718 b = sub_b->a; 3719 c = sub_c->v; 3720 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3721 3722 ii = sub_b->i; 3723 idx = sub_b->j; 3724 for (i=0; i<n; i++) { 3725 q = ii[i+1] - ii[i]; 3726 while (q-->0) { 3727 c_q = c + m*(*idx); 3728 a_q = a + m*i; 3729 PetscKernelAXPY(c_q,*b,a_q,m); 3730 idx++; 3731 b++; 3732 } 3733 } 3734 PetscFunctionReturn(0); 3735 } 3736 3737 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3738 { 3739 PetscErrorCode ierr; 3740 PetscInt m=A->rmap->n,n=B->cmap->n; 3741 Mat Cmat; 3742 3743 PetscFunctionBegin; 3744 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); 3745 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3746 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3747 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3748 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3749 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3750 3751 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3752 3753 *C = Cmat; 3754 PetscFunctionReturn(0); 3755 } 3756 3757 /* ----------------------------------------------------------------*/ 3758 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3759 { 3760 PetscErrorCode ierr; 3761 3762 PetscFunctionBegin; 3763 if (scall == MAT_INITIAL_MATRIX) { 3764 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3765 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3766 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3767 } 3768 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3769 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3770 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3771 PetscFunctionReturn(0); 3772 } 3773 3774 3775 /*MC 3776 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3777 based on compressed sparse row format. 3778 3779 Options Database Keys: 3780 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3781 3782 Level: beginner 3783 3784 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3785 M*/ 3786 3787 /*MC 3788 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3789 3790 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3791 and MATMPIAIJ otherwise. As a result, for single process communicators, 3792 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3793 for communicators controlling multiple processes. It is recommended that you call both of 3794 the above preallocation routines for simplicity. 3795 3796 Options Database Keys: 3797 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3798 3799 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3800 enough exist. 3801 3802 Level: beginner 3803 3804 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3805 M*/ 3806 3807 /*MC 3808 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3809 3810 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3811 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3812 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3813 for communicators controlling multiple processes. It is recommended that you call both of 3814 the above preallocation routines for simplicity. 3815 3816 Options Database Keys: 3817 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3818 3819 Level: beginner 3820 3821 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3822 M*/ 3823 3824 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3825 #if defined(PETSC_HAVE_ELEMENTAL) 3826 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3827 #endif 3828 #if defined(PETSC_HAVE_HYPRE) 3829 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3830 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3831 #endif 3832 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3833 3834 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3835 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3836 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3837 #endif 3838 3839 3840 /*@C 3841 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 3842 3843 Not Collective 3844 3845 Input Parameter: 3846 . mat - a MATSEQAIJ matrix 3847 3848 Output Parameter: 3849 . array - pointer to the data 3850 3851 Level: intermediate 3852 3853 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3854 @*/ 3855 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3856 { 3857 PetscErrorCode ierr; 3858 3859 PetscFunctionBegin; 3860 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3861 PetscFunctionReturn(0); 3862 } 3863 3864 /*@C 3865 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3866 3867 Not Collective 3868 3869 Input Parameter: 3870 . mat - a MATSEQAIJ matrix 3871 3872 Output Parameter: 3873 . nz - the maximum number of nonzeros in any row 3874 3875 Level: intermediate 3876 3877 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3878 @*/ 3879 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3880 { 3881 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3882 3883 PetscFunctionBegin; 3884 *nz = aij->rmax; 3885 PetscFunctionReturn(0); 3886 } 3887 3888 /*@C 3889 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3890 3891 Not Collective 3892 3893 Input Parameters: 3894 . mat - a MATSEQAIJ matrix 3895 . array - pointer to the data 3896 3897 Level: intermediate 3898 3899 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3900 @*/ 3901 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3902 { 3903 PetscErrorCode ierr; 3904 3905 PetscFunctionBegin; 3906 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3907 PetscFunctionReturn(0); 3908 } 3909 3910 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3911 { 3912 Mat_SeqAIJ *b; 3913 PetscErrorCode ierr; 3914 PetscMPIInt size; 3915 3916 PetscFunctionBegin; 3917 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3918 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3919 3920 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3921 3922 B->data = (void*)b; 3923 3924 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3925 3926 b->row = 0; 3927 b->col = 0; 3928 b->icol = 0; 3929 b->reallocs = 0; 3930 b->ignorezeroentries = PETSC_FALSE; 3931 b->roworiented = PETSC_TRUE; 3932 b->nonew = 0; 3933 b->diag = 0; 3934 b->solve_work = 0; 3935 B->spptr = 0; 3936 b->saved_values = 0; 3937 b->idiag = 0; 3938 b->mdiag = 0; 3939 b->ssor_work = 0; 3940 b->omega = 1.0; 3941 b->fshift = 0.0; 3942 b->idiagvalid = PETSC_FALSE; 3943 b->ibdiagvalid = PETSC_FALSE; 3944 b->keepnonzeropattern = PETSC_FALSE; 3945 3946 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3947 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3948 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3949 3950 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3951 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3952 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3953 #endif 3954 3955 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3956 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3957 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3958 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3959 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3960 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3961 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3962 #if defined(PETSC_HAVE_ELEMENTAL) 3963 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 3964 #endif 3965 #if defined(PETSC_HAVE_HYPRE) 3966 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 3967 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 3968 #endif 3969 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 3970 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3971 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3972 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3973 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3974 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3975 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3976 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3977 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3978 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3979 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3980 PetscFunctionReturn(0); 3981 } 3982 3983 /* 3984 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3985 */ 3986 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3987 { 3988 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3989 PetscErrorCode ierr; 3990 PetscInt i,m = A->rmap->n; 3991 3992 PetscFunctionBegin; 3993 c = (Mat_SeqAIJ*)C->data; 3994 3995 C->factortype = A->factortype; 3996 c->row = 0; 3997 c->col = 0; 3998 c->icol = 0; 3999 c->reallocs = 0; 4000 4001 C->assembled = PETSC_TRUE; 4002 4003 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4004 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4005 4006 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 4007 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4008 for (i=0; i<m; i++) { 4009 c->imax[i] = a->imax[i]; 4010 c->ilen[i] = a->ilen[i]; 4011 } 4012 4013 /* allocate the matrix space */ 4014 if (mallocmatspace) { 4015 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 4016 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4017 4018 c->singlemalloc = PETSC_TRUE; 4019 4020 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4021 if (m > 0) { 4022 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4023 if (cpvalues == MAT_COPY_VALUES) { 4024 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4025 } else { 4026 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4027 } 4028 } 4029 } 4030 4031 c->ignorezeroentries = a->ignorezeroentries; 4032 c->roworiented = a->roworiented; 4033 c->nonew = a->nonew; 4034 if (a->diag) { 4035 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4036 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4037 for (i=0; i<m; i++) { 4038 c->diag[i] = a->diag[i]; 4039 } 4040 } else c->diag = 0; 4041 4042 c->solve_work = 0; 4043 c->saved_values = 0; 4044 c->idiag = 0; 4045 c->ssor_work = 0; 4046 c->keepnonzeropattern = a->keepnonzeropattern; 4047 c->free_a = PETSC_TRUE; 4048 c->free_ij = PETSC_TRUE; 4049 4050 c->rmax = a->rmax; 4051 c->nz = a->nz; 4052 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4053 C->preallocated = PETSC_TRUE; 4054 4055 c->compressedrow.use = a->compressedrow.use; 4056 c->compressedrow.nrows = a->compressedrow.nrows; 4057 if (a->compressedrow.use) { 4058 i = a->compressedrow.nrows; 4059 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4060 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4061 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4062 } else { 4063 c->compressedrow.use = PETSC_FALSE; 4064 c->compressedrow.i = NULL; 4065 c->compressedrow.rindex = NULL; 4066 } 4067 c->nonzerorowcnt = a->nonzerorowcnt; 4068 C->nonzerostate = A->nonzerostate; 4069 4070 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4071 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4072 PetscFunctionReturn(0); 4073 } 4074 4075 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4076 { 4077 PetscErrorCode ierr; 4078 4079 PetscFunctionBegin; 4080 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4081 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4082 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4083 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4084 } 4085 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4086 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4087 PetscFunctionReturn(0); 4088 } 4089 4090 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4091 { 4092 Mat_SeqAIJ *a; 4093 PetscErrorCode ierr; 4094 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4095 int fd; 4096 PetscMPIInt size; 4097 MPI_Comm comm; 4098 PetscInt bs = newMat->rmap->bs; 4099 4100 PetscFunctionBegin; 4101 /* force binary viewer to load .info file if it has not yet done so */ 4102 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4103 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4104 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4105 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4106 4107 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4108 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4109 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4110 if (bs < 0) bs = 1; 4111 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4112 4113 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4114 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4115 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4116 M = header[1]; N = header[2]; nz = header[3]; 4117 4118 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4119 4120 /* read in row lengths */ 4121 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4122 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4123 4124 /* check if sum of rowlengths is same as nz */ 4125 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4126 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); 4127 4128 /* set global size if not set already*/ 4129 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4130 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4131 } else { 4132 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4133 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4134 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4135 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4136 } 4137 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); 4138 } 4139 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4140 a = (Mat_SeqAIJ*)newMat->data; 4141 4142 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4143 4144 /* read in nonzero values */ 4145 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4146 4147 /* set matrix "i" values */ 4148 a->i[0] = 0; 4149 for (i=1; i<= M; i++) { 4150 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4151 a->ilen[i-1] = rowlengths[i-1]; 4152 } 4153 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4154 4155 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4156 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4157 PetscFunctionReturn(0); 4158 } 4159 4160 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4161 { 4162 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4163 PetscErrorCode ierr; 4164 #if defined(PETSC_USE_COMPLEX) 4165 PetscInt k; 4166 #endif 4167 4168 PetscFunctionBegin; 4169 /* If the matrix dimensions are not equal,or no of nonzeros */ 4170 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4171 *flg = PETSC_FALSE; 4172 PetscFunctionReturn(0); 4173 } 4174 4175 /* if the a->i are the same */ 4176 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4177 if (!*flg) PetscFunctionReturn(0); 4178 4179 /* if a->j are the same */ 4180 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4181 if (!*flg) PetscFunctionReturn(0); 4182 4183 /* if a->a are the same */ 4184 #if defined(PETSC_USE_COMPLEX) 4185 for (k=0; k<a->nz; k++) { 4186 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4187 *flg = PETSC_FALSE; 4188 PetscFunctionReturn(0); 4189 } 4190 } 4191 #else 4192 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4193 #endif 4194 PetscFunctionReturn(0); 4195 } 4196 4197 /*@ 4198 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4199 provided by the user. 4200 4201 Collective on MPI_Comm 4202 4203 Input Parameters: 4204 + comm - must be an MPI communicator of size 1 4205 . m - number of rows 4206 . n - number of columns 4207 . i - row indices 4208 . j - column indices 4209 - a - matrix values 4210 4211 Output Parameter: 4212 . mat - the matrix 4213 4214 Level: intermediate 4215 4216 Notes: 4217 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4218 once the matrix is destroyed and not before 4219 4220 You cannot set new nonzero locations into this matrix, that will generate an error. 4221 4222 The i and j indices are 0 based 4223 4224 The format which is used for the sparse matrix input, is equivalent to a 4225 row-major ordering.. i.e for the following matrix, the input data expected is 4226 as shown 4227 4228 $ 1 0 0 4229 $ 2 0 3 4230 $ 4 5 6 4231 $ 4232 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4233 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4234 $ v = {1,2,3,4,5,6} [size = 6] 4235 4236 4237 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4238 4239 @*/ 4240 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4241 { 4242 PetscErrorCode ierr; 4243 PetscInt ii; 4244 Mat_SeqAIJ *aij; 4245 #if defined(PETSC_USE_DEBUG) 4246 PetscInt jj; 4247 #endif 4248 4249 PetscFunctionBegin; 4250 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4251 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4252 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4253 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4254 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4255 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4256 aij = (Mat_SeqAIJ*)(*mat)->data; 4257 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4258 4259 aij->i = i; 4260 aij->j = j; 4261 aij->a = a; 4262 aij->singlemalloc = PETSC_FALSE; 4263 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4264 aij->free_a = PETSC_FALSE; 4265 aij->free_ij = PETSC_FALSE; 4266 4267 for (ii=0; ii<m; ii++) { 4268 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4269 #if defined(PETSC_USE_DEBUG) 4270 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]); 4271 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 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 not sorted",jj-i[ii],j[jj],ii); 4273 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); 4274 } 4275 #endif 4276 } 4277 #if defined(PETSC_USE_DEBUG) 4278 for (ii=0; ii<aij->i[m]; ii++) { 4279 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4280 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]); 4281 } 4282 #endif 4283 4284 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4285 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4286 PetscFunctionReturn(0); 4287 } 4288 /*@C 4289 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4290 provided by the user. 4291 4292 Collective on MPI_Comm 4293 4294 Input Parameters: 4295 + comm - must be an MPI communicator of size 1 4296 . m - number of rows 4297 . n - number of columns 4298 . i - row indices 4299 . j - column indices 4300 . a - matrix values 4301 . nz - number of nonzeros 4302 - idx - 0 or 1 based 4303 4304 Output Parameter: 4305 . mat - the matrix 4306 4307 Level: intermediate 4308 4309 Notes: 4310 The i and j indices are 0 based 4311 4312 The format which is used for the sparse matrix input, is equivalent to a 4313 row-major ordering.. i.e for the following matrix, the input data expected is 4314 as shown: 4315 4316 1 0 0 4317 2 0 3 4318 4 5 6 4319 4320 i = {0,1,1,2,2,2} 4321 j = {0,0,2,0,1,2} 4322 v = {1,2,3,4,5,6} 4323 4324 4325 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4326 4327 @*/ 4328 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4329 { 4330 PetscErrorCode ierr; 4331 PetscInt ii, *nnz, one = 1,row,col; 4332 4333 4334 PetscFunctionBegin; 4335 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4336 for (ii = 0; ii < nz; ii++) { 4337 nnz[i[ii] - !!idx] += 1; 4338 } 4339 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4340 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4341 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4342 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4343 for (ii = 0; ii < nz; ii++) { 4344 if (idx) { 4345 row = i[ii] - 1; 4346 col = j[ii] - 1; 4347 } else { 4348 row = i[ii]; 4349 col = j[ii]; 4350 } 4351 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4352 } 4353 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4354 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4355 ierr = PetscFree(nnz);CHKERRQ(ierr); 4356 PetscFunctionReturn(0); 4357 } 4358 4359 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4360 { 4361 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4362 PetscErrorCode ierr; 4363 4364 PetscFunctionBegin; 4365 a->idiagvalid = PETSC_FALSE; 4366 a->ibdiagvalid = PETSC_FALSE; 4367 4368 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4369 PetscFunctionReturn(0); 4370 } 4371 4372 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4373 { 4374 PetscErrorCode ierr; 4375 PetscMPIInt size; 4376 4377 PetscFunctionBegin; 4378 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4379 if (size == 1 && scall == MAT_REUSE_MATRIX) { 4380 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4381 } else { 4382 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4383 } 4384 PetscFunctionReturn(0); 4385 } 4386 4387 /* 4388 Permute A into C's *local* index space using rowemb,colemb. 4389 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4390 of [0,m), colemb is in [0,n). 4391 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4392 */ 4393 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4394 { 4395 /* If making this function public, change the error returned in this function away from _PLIB. */ 4396 PetscErrorCode ierr; 4397 Mat_SeqAIJ *Baij; 4398 PetscBool seqaij; 4399 PetscInt m,n,*nz,i,j,count; 4400 PetscScalar v; 4401 const PetscInt *rowindices,*colindices; 4402 4403 PetscFunctionBegin; 4404 if (!B) PetscFunctionReturn(0); 4405 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4406 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4407 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4408 if (rowemb) { 4409 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4410 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); 4411 } else { 4412 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4413 } 4414 if (colemb) { 4415 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4416 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); 4417 } else { 4418 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4419 } 4420 4421 Baij = (Mat_SeqAIJ*)(B->data); 4422 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4423 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4424 for (i=0; i<B->rmap->n; i++) { 4425 nz[i] = Baij->i[i+1] - Baij->i[i]; 4426 } 4427 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4428 ierr = PetscFree(nz);CHKERRQ(ierr); 4429 } 4430 if (pattern == SUBSET_NONZERO_PATTERN) { 4431 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4432 } 4433 count = 0; 4434 rowindices = NULL; 4435 colindices = NULL; 4436 if (rowemb) { 4437 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4438 } 4439 if (colemb) { 4440 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4441 } 4442 for (i=0; i<B->rmap->n; i++) { 4443 PetscInt row; 4444 row = i; 4445 if (rowindices) row = rowindices[i]; 4446 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4447 PetscInt col; 4448 col = Baij->j[count]; 4449 if (colindices) col = colindices[col]; 4450 v = Baij->a[count]; 4451 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4452 ++count; 4453 } 4454 } 4455 /* FIXME: set C's nonzerostate correctly. */ 4456 /* Assembly for C is necessary. */ 4457 C->preallocated = PETSC_TRUE; 4458 C->assembled = PETSC_TRUE; 4459 C->was_assembled = PETSC_FALSE; 4460 PetscFunctionReturn(0); 4461 } 4462 4463 4464 /* 4465 Special version for direct calls from Fortran 4466 */ 4467 #include <petsc/private/fortranimpl.h> 4468 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4469 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4470 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4471 #define matsetvaluesseqaij_ matsetvaluesseqaij 4472 #endif 4473 4474 /* Change these macros so can be used in void function */ 4475 #undef CHKERRQ 4476 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4477 #undef SETERRQ2 4478 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4479 #undef SETERRQ3 4480 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4481 4482 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) 4483 { 4484 Mat A = *AA; 4485 PetscInt m = *mm, n = *nn; 4486 InsertMode is = *isis; 4487 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4488 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4489 PetscInt *imax,*ai,*ailen; 4490 PetscErrorCode ierr; 4491 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4492 MatScalar *ap,value,*aa; 4493 PetscBool ignorezeroentries = a->ignorezeroentries; 4494 PetscBool roworiented = a->roworiented; 4495 4496 PetscFunctionBegin; 4497 MatCheckPreallocated(A,1); 4498 imax = a->imax; 4499 ai = a->i; 4500 ailen = a->ilen; 4501 aj = a->j; 4502 aa = a->a; 4503 4504 for (k=0; k<m; k++) { /* loop over added rows */ 4505 row = im[k]; 4506 if (row < 0) continue; 4507 #if defined(PETSC_USE_DEBUG) 4508 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4509 #endif 4510 rp = aj + ai[row]; ap = aa + ai[row]; 4511 rmax = imax[row]; nrow = ailen[row]; 4512 low = 0; 4513 high = nrow; 4514 for (l=0; l<n; l++) { /* loop over added columns */ 4515 if (in[l] < 0) continue; 4516 #if defined(PETSC_USE_DEBUG) 4517 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4518 #endif 4519 col = in[l]; 4520 if (roworiented) value = v[l + k*n]; 4521 else value = v[k + l*m]; 4522 4523 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4524 4525 if (col <= lastcol) low = 0; 4526 else high = nrow; 4527 lastcol = col; 4528 while (high-low > 5) { 4529 t = (low+high)/2; 4530 if (rp[t] > col) high = t; 4531 else low = t; 4532 } 4533 for (i=low; i<high; i++) { 4534 if (rp[i] > col) break; 4535 if (rp[i] == col) { 4536 if (is == ADD_VALUES) ap[i] += value; 4537 else ap[i] = value; 4538 goto noinsert; 4539 } 4540 } 4541 if (value == 0.0 && ignorezeroentries) goto noinsert; 4542 if (nonew == 1) goto noinsert; 4543 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4544 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4545 N = nrow++ - 1; a->nz++; high++; 4546 /* shift up all the later entries in this row */ 4547 for (ii=N; ii>=i; ii--) { 4548 rp[ii+1] = rp[ii]; 4549 ap[ii+1] = ap[ii]; 4550 } 4551 rp[i] = col; 4552 ap[i] = value; 4553 A->nonzerostate++; 4554 noinsert:; 4555 low = i + 1; 4556 } 4557 ailen[row] = nrow; 4558 } 4559 PetscFunctionReturnVoid(); 4560 } 4561 4562