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