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