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