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 = PetscMalloc2(nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3558 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*sizeof(PetscInt));CHKERRQ(ierr); 3559 } else { 3560 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3561 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3562 } 3563 b->i[0] = 0; 3564 for (i=1; i<B->rmap->n+1; i++) { 3565 b->i[i] = b->i[i-1] + b->imax[i-1]; 3566 } 3567 if (B->structure_only) { 3568 b->singlemalloc = PETSC_FALSE; 3569 b->free_a = PETSC_FALSE; 3570 } else { 3571 b->singlemalloc = PETSC_TRUE; 3572 b->free_a = PETSC_TRUE; 3573 } 3574 b->free_ij = PETSC_TRUE; 3575 } else { 3576 b->free_a = PETSC_FALSE; 3577 b->free_ij = PETSC_FALSE; 3578 } 3579 3580 b->nz = 0; 3581 b->maxnz = nz; 3582 B->info.nz_unneeded = (double)b->maxnz; 3583 if (realalloc) { 3584 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3585 } 3586 B->was_assembled = PETSC_FALSE; 3587 B->assembled = PETSC_FALSE; 3588 PetscFunctionReturn(0); 3589 } 3590 3591 /*@ 3592 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3593 3594 Input Parameters: 3595 + B - the matrix 3596 . i - the indices into j for the start of each row (starts with zero) 3597 . j - the column indices for each row (starts with zero) these must be sorted for each row 3598 - v - optional values in the matrix 3599 3600 Level: developer 3601 3602 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3603 3604 .keywords: matrix, aij, compressed row, sparse, sequential 3605 3606 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3607 @*/ 3608 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3609 { 3610 PetscErrorCode ierr; 3611 3612 PetscFunctionBegin; 3613 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3614 PetscValidType(B,1); 3615 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3616 PetscFunctionReturn(0); 3617 } 3618 3619 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3620 { 3621 PetscInt i; 3622 PetscInt m,n; 3623 PetscInt nz; 3624 PetscInt *nnz, nz_max = 0; 3625 PetscScalar *values; 3626 PetscErrorCode ierr; 3627 3628 PetscFunctionBegin; 3629 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3630 3631 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3632 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3633 3634 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3635 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3636 for (i = 0; i < m; i++) { 3637 nz = Ii[i+1]- Ii[i]; 3638 nz_max = PetscMax(nz_max, nz); 3639 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3640 nnz[i] = nz; 3641 } 3642 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3643 ierr = PetscFree(nnz);CHKERRQ(ierr); 3644 3645 if (v) { 3646 values = (PetscScalar*) v; 3647 } else { 3648 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3649 } 3650 3651 for (i = 0; i < m; i++) { 3652 nz = Ii[i+1] - Ii[i]; 3653 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3654 } 3655 3656 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3657 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3658 3659 if (!v) { 3660 ierr = PetscFree(values);CHKERRQ(ierr); 3661 } 3662 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3663 PetscFunctionReturn(0); 3664 } 3665 3666 #include <../src/mat/impls/dense/seq/dense.h> 3667 #include <petsc/private/kernels/petscaxpy.h> 3668 3669 /* 3670 Computes (B'*A')' since computing B*A directly is untenable 3671 3672 n p p 3673 ( ) ( ) ( ) 3674 m ( A ) * n ( B ) = m ( C ) 3675 ( ) ( ) ( ) 3676 3677 */ 3678 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3679 { 3680 PetscErrorCode ierr; 3681 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3682 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3683 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3684 PetscInt i,n,m,q,p; 3685 const PetscInt *ii,*idx; 3686 const PetscScalar *b,*a,*a_q; 3687 PetscScalar *c,*c_q; 3688 3689 PetscFunctionBegin; 3690 m = A->rmap->n; 3691 n = A->cmap->n; 3692 p = B->cmap->n; 3693 a = sub_a->v; 3694 b = sub_b->a; 3695 c = sub_c->v; 3696 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3697 3698 ii = sub_b->i; 3699 idx = sub_b->j; 3700 for (i=0; i<n; i++) { 3701 q = ii[i+1] - ii[i]; 3702 while (q-->0) { 3703 c_q = c + m*(*idx); 3704 a_q = a + m*i; 3705 PetscKernelAXPY(c_q,*b,a_q,m); 3706 idx++; 3707 b++; 3708 } 3709 } 3710 PetscFunctionReturn(0); 3711 } 3712 3713 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3714 { 3715 PetscErrorCode ierr; 3716 PetscInt m=A->rmap->n,n=B->cmap->n; 3717 Mat Cmat; 3718 3719 PetscFunctionBegin; 3720 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); 3721 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3722 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3723 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3724 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3725 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3726 3727 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3728 3729 *C = Cmat; 3730 PetscFunctionReturn(0); 3731 } 3732 3733 /* ----------------------------------------------------------------*/ 3734 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3735 { 3736 PetscErrorCode ierr; 3737 3738 PetscFunctionBegin; 3739 if (scall == MAT_INITIAL_MATRIX) { 3740 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3741 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3742 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3743 } 3744 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3745 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3746 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3747 PetscFunctionReturn(0); 3748 } 3749 3750 3751 /*MC 3752 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3753 based on compressed sparse row format. 3754 3755 Options Database Keys: 3756 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3757 3758 Level: beginner 3759 3760 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3761 M*/ 3762 3763 /*MC 3764 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3765 3766 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3767 and MATMPIAIJ otherwise. As a result, for single process communicators, 3768 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3769 for communicators controlling multiple processes. It is recommended that you call both of 3770 the above preallocation routines for simplicity. 3771 3772 Options Database Keys: 3773 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3774 3775 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3776 enough exist. 3777 3778 Level: beginner 3779 3780 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3781 M*/ 3782 3783 /*MC 3784 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3785 3786 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3787 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3788 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3789 for communicators controlling multiple processes. It is recommended that you call both of 3790 the above preallocation routines for simplicity. 3791 3792 Options Database Keys: 3793 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3794 3795 Level: beginner 3796 3797 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3798 M*/ 3799 3800 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3801 #if defined(PETSC_HAVE_ELEMENTAL) 3802 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3803 #endif 3804 #if defined(PETSC_HAVE_HYPRE) 3805 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3806 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3807 #endif 3808 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3809 3810 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3811 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3812 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3813 #endif 3814 3815 3816 /*@C 3817 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 3818 3819 Not Collective 3820 3821 Input Parameter: 3822 . mat - a MATSEQAIJ matrix 3823 3824 Output Parameter: 3825 . array - pointer to the data 3826 3827 Level: intermediate 3828 3829 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3830 @*/ 3831 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3832 { 3833 PetscErrorCode ierr; 3834 3835 PetscFunctionBegin; 3836 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3837 PetscFunctionReturn(0); 3838 } 3839 3840 /*@C 3841 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3842 3843 Not Collective 3844 3845 Input Parameter: 3846 . mat - a MATSEQAIJ matrix 3847 3848 Output Parameter: 3849 . nz - the maximum number of nonzeros in any row 3850 3851 Level: intermediate 3852 3853 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3854 @*/ 3855 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3856 { 3857 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3858 3859 PetscFunctionBegin; 3860 *nz = aij->rmax; 3861 PetscFunctionReturn(0); 3862 } 3863 3864 /*@C 3865 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3866 3867 Not Collective 3868 3869 Input Parameters: 3870 . mat - a MATSEQAIJ matrix 3871 . array - pointer to the data 3872 3873 Level: intermediate 3874 3875 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3876 @*/ 3877 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3878 { 3879 PetscErrorCode ierr; 3880 3881 PetscFunctionBegin; 3882 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3883 PetscFunctionReturn(0); 3884 } 3885 3886 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3887 { 3888 Mat_SeqAIJ *b; 3889 PetscErrorCode ierr; 3890 PetscMPIInt size; 3891 3892 PetscFunctionBegin; 3893 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3894 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3895 3896 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3897 3898 B->data = (void*)b; 3899 3900 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3901 3902 b->row = 0; 3903 b->col = 0; 3904 b->icol = 0; 3905 b->reallocs = 0; 3906 b->ignorezeroentries = PETSC_FALSE; 3907 b->roworiented = PETSC_TRUE; 3908 b->nonew = 0; 3909 b->diag = 0; 3910 b->solve_work = 0; 3911 B->spptr = 0; 3912 b->saved_values = 0; 3913 b->idiag = 0; 3914 b->mdiag = 0; 3915 b->ssor_work = 0; 3916 b->omega = 1.0; 3917 b->fshift = 0.0; 3918 b->idiagvalid = PETSC_FALSE; 3919 b->ibdiagvalid = PETSC_FALSE; 3920 b->keepnonzeropattern = PETSC_FALSE; 3921 3922 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3923 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3924 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3925 3926 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3927 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3928 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3929 #endif 3930 3931 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3932 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3933 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3934 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3935 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3936 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3937 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3938 #if defined(PETSC_HAVE_ELEMENTAL) 3939 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 3940 #endif 3941 #if defined(PETSC_HAVE_HYPRE) 3942 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 3943 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 3944 #endif 3945 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 3946 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3947 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3948 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3949 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3950 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3951 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3952 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3953 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3954 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3955 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3956 PetscFunctionReturn(0); 3957 } 3958 3959 /* 3960 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3961 */ 3962 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3963 { 3964 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3965 PetscErrorCode ierr; 3966 PetscInt i,m = A->rmap->n; 3967 3968 PetscFunctionBegin; 3969 c = (Mat_SeqAIJ*)C->data; 3970 3971 C->factortype = A->factortype; 3972 c->row = 0; 3973 c->col = 0; 3974 c->icol = 0; 3975 c->reallocs = 0; 3976 3977 C->assembled = PETSC_TRUE; 3978 3979 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3980 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3981 3982 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 3983 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3984 for (i=0; i<m; i++) { 3985 c->imax[i] = a->imax[i]; 3986 c->ilen[i] = a->ilen[i]; 3987 } 3988 3989 /* allocate the matrix space */ 3990 if (mallocmatspace) { 3991 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 3992 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3993 3994 c->singlemalloc = PETSC_TRUE; 3995 3996 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3997 if (m > 0) { 3998 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3999 if (cpvalues == MAT_COPY_VALUES) { 4000 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4001 } else { 4002 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4003 } 4004 } 4005 } 4006 4007 c->ignorezeroentries = a->ignorezeroentries; 4008 c->roworiented = a->roworiented; 4009 c->nonew = a->nonew; 4010 if (a->diag) { 4011 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 4012 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4013 for (i=0; i<m; i++) { 4014 c->diag[i] = a->diag[i]; 4015 } 4016 } else c->diag = 0; 4017 4018 c->solve_work = 0; 4019 c->saved_values = 0; 4020 c->idiag = 0; 4021 c->ssor_work = 0; 4022 c->keepnonzeropattern = a->keepnonzeropattern; 4023 c->free_a = PETSC_TRUE; 4024 c->free_ij = PETSC_TRUE; 4025 4026 c->rmax = a->rmax; 4027 c->nz = a->nz; 4028 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4029 C->preallocated = PETSC_TRUE; 4030 4031 c->compressedrow.use = a->compressedrow.use; 4032 c->compressedrow.nrows = a->compressedrow.nrows; 4033 if (a->compressedrow.use) { 4034 i = a->compressedrow.nrows; 4035 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4036 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4037 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4038 } else { 4039 c->compressedrow.use = PETSC_FALSE; 4040 c->compressedrow.i = NULL; 4041 c->compressedrow.rindex = NULL; 4042 } 4043 c->nonzerorowcnt = a->nonzerorowcnt; 4044 C->nonzerostate = A->nonzerostate; 4045 4046 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4047 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4048 PetscFunctionReturn(0); 4049 } 4050 4051 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4052 { 4053 PetscErrorCode ierr; 4054 4055 PetscFunctionBegin; 4056 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4057 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4058 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4059 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4060 } 4061 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4062 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4063 PetscFunctionReturn(0); 4064 } 4065 4066 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4067 { 4068 Mat_SeqAIJ *a; 4069 PetscErrorCode ierr; 4070 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4071 int fd; 4072 PetscMPIInt size; 4073 MPI_Comm comm; 4074 PetscInt bs = newMat->rmap->bs; 4075 4076 PetscFunctionBegin; 4077 /* force binary viewer to load .info file if it has not yet done so */ 4078 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4079 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4080 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4081 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4082 4083 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4084 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4085 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4086 if (bs < 0) bs = 1; 4087 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4088 4089 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4090 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4091 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4092 M = header[1]; N = header[2]; nz = header[3]; 4093 4094 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4095 4096 /* read in row lengths */ 4097 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4098 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4099 4100 /* check if sum of rowlengths is same as nz */ 4101 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4102 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); 4103 4104 /* set global size if not set already*/ 4105 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4106 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4107 } else { 4108 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4109 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4110 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4111 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4112 } 4113 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); 4114 } 4115 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4116 a = (Mat_SeqAIJ*)newMat->data; 4117 4118 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4119 4120 /* read in nonzero values */ 4121 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4122 4123 /* set matrix "i" values */ 4124 a->i[0] = 0; 4125 for (i=1; i<= M; i++) { 4126 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4127 a->ilen[i-1] = rowlengths[i-1]; 4128 } 4129 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4130 4131 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4132 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4133 PetscFunctionReturn(0); 4134 } 4135 4136 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4137 { 4138 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4139 PetscErrorCode ierr; 4140 #if defined(PETSC_USE_COMPLEX) 4141 PetscInt k; 4142 #endif 4143 4144 PetscFunctionBegin; 4145 /* If the matrix dimensions are not equal,or no of nonzeros */ 4146 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4147 *flg = PETSC_FALSE; 4148 PetscFunctionReturn(0); 4149 } 4150 4151 /* if the a->i are the same */ 4152 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4153 if (!*flg) PetscFunctionReturn(0); 4154 4155 /* if a->j are the same */ 4156 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4157 if (!*flg) PetscFunctionReturn(0); 4158 4159 /* if a->a are the same */ 4160 #if defined(PETSC_USE_COMPLEX) 4161 for (k=0; k<a->nz; k++) { 4162 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4163 *flg = PETSC_FALSE; 4164 PetscFunctionReturn(0); 4165 } 4166 } 4167 #else 4168 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4169 #endif 4170 PetscFunctionReturn(0); 4171 } 4172 4173 /*@ 4174 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4175 provided by the user. 4176 4177 Collective on MPI_Comm 4178 4179 Input Parameters: 4180 + comm - must be an MPI communicator of size 1 4181 . m - number of rows 4182 . n - number of columns 4183 . i - row indices 4184 . j - column indices 4185 - a - matrix values 4186 4187 Output Parameter: 4188 . mat - the matrix 4189 4190 Level: intermediate 4191 4192 Notes: 4193 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4194 once the matrix is destroyed and not before 4195 4196 You cannot set new nonzero locations into this matrix, that will generate an error. 4197 4198 The i and j indices are 0 based 4199 4200 The format which is used for the sparse matrix input, is equivalent to a 4201 row-major ordering.. i.e for the following matrix, the input data expected is 4202 as shown 4203 4204 $ 1 0 0 4205 $ 2 0 3 4206 $ 4 5 6 4207 $ 4208 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4209 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4210 $ v = {1,2,3,4,5,6} [size = 6] 4211 4212 4213 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4214 4215 @*/ 4216 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4217 { 4218 PetscErrorCode ierr; 4219 PetscInt ii; 4220 Mat_SeqAIJ *aij; 4221 #if defined(PETSC_USE_DEBUG) 4222 PetscInt jj; 4223 #endif 4224 4225 PetscFunctionBegin; 4226 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4227 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4228 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4229 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4230 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4231 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4232 aij = (Mat_SeqAIJ*)(*mat)->data; 4233 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4234 4235 aij->i = i; 4236 aij->j = j; 4237 aij->a = a; 4238 aij->singlemalloc = PETSC_FALSE; 4239 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4240 aij->free_a = PETSC_FALSE; 4241 aij->free_ij = PETSC_FALSE; 4242 4243 for (ii=0; ii<m; ii++) { 4244 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4245 #if defined(PETSC_USE_DEBUG) 4246 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]); 4247 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4248 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); 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 identical to previous entry",jj-i[ii],j[jj],ii); 4250 } 4251 #endif 4252 } 4253 #if defined(PETSC_USE_DEBUG) 4254 for (ii=0; ii<aij->i[m]; ii++) { 4255 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4256 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]); 4257 } 4258 #endif 4259 4260 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4261 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4262 PetscFunctionReturn(0); 4263 } 4264 /*@C 4265 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4266 provided by the user. 4267 4268 Collective on MPI_Comm 4269 4270 Input Parameters: 4271 + comm - must be an MPI communicator of size 1 4272 . m - number of rows 4273 . n - number of columns 4274 . i - row indices 4275 . j - column indices 4276 . a - matrix values 4277 . nz - number of nonzeros 4278 - idx - 0 or 1 based 4279 4280 Output Parameter: 4281 . mat - the matrix 4282 4283 Level: intermediate 4284 4285 Notes: 4286 The i and j indices are 0 based 4287 4288 The format which is used for the sparse matrix input, is equivalent to a 4289 row-major ordering.. i.e for the following matrix, the input data expected is 4290 as shown: 4291 4292 1 0 0 4293 2 0 3 4294 4 5 6 4295 4296 i = {0,1,1,2,2,2} 4297 j = {0,0,2,0,1,2} 4298 v = {1,2,3,4,5,6} 4299 4300 4301 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4302 4303 @*/ 4304 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4305 { 4306 PetscErrorCode ierr; 4307 PetscInt ii, *nnz, one = 1,row,col; 4308 4309 4310 PetscFunctionBegin; 4311 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4312 for (ii = 0; ii < nz; ii++) { 4313 nnz[i[ii] - !!idx] += 1; 4314 } 4315 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4316 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4317 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4318 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4319 for (ii = 0; ii < nz; ii++) { 4320 if (idx) { 4321 row = i[ii] - 1; 4322 col = j[ii] - 1; 4323 } else { 4324 row = i[ii]; 4325 col = j[ii]; 4326 } 4327 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4328 } 4329 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4330 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4331 ierr = PetscFree(nnz);CHKERRQ(ierr); 4332 PetscFunctionReturn(0); 4333 } 4334 4335 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4336 { 4337 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4338 PetscErrorCode ierr; 4339 4340 PetscFunctionBegin; 4341 a->idiagvalid = PETSC_FALSE; 4342 a->ibdiagvalid = PETSC_FALSE; 4343 4344 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4345 PetscFunctionReturn(0); 4346 } 4347 4348 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4349 { 4350 PetscErrorCode ierr; 4351 PetscMPIInt size; 4352 4353 PetscFunctionBegin; 4354 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4355 if (size == 1 && scall == MAT_REUSE_MATRIX) { 4356 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4357 } else { 4358 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4359 } 4360 PetscFunctionReturn(0); 4361 } 4362 4363 /* 4364 Permute A into C's *local* index space using rowemb,colemb. 4365 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4366 of [0,m), colemb is in [0,n). 4367 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4368 */ 4369 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4370 { 4371 /* If making this function public, change the error returned in this function away from _PLIB. */ 4372 PetscErrorCode ierr; 4373 Mat_SeqAIJ *Baij; 4374 PetscBool seqaij; 4375 PetscInt m,n,*nz,i,j,count; 4376 PetscScalar v; 4377 const PetscInt *rowindices,*colindices; 4378 4379 PetscFunctionBegin; 4380 if (!B) PetscFunctionReturn(0); 4381 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4382 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4383 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4384 if (rowemb) { 4385 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4386 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); 4387 } else { 4388 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4389 } 4390 if (colemb) { 4391 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4392 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); 4393 } else { 4394 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4395 } 4396 4397 Baij = (Mat_SeqAIJ*)(B->data); 4398 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4399 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4400 for (i=0; i<B->rmap->n; i++) { 4401 nz[i] = Baij->i[i+1] - Baij->i[i]; 4402 } 4403 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4404 ierr = PetscFree(nz);CHKERRQ(ierr); 4405 } 4406 if (pattern == SUBSET_NONZERO_PATTERN) { 4407 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4408 } 4409 count = 0; 4410 rowindices = NULL; 4411 colindices = NULL; 4412 if (rowemb) { 4413 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4414 } 4415 if (colemb) { 4416 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4417 } 4418 for (i=0; i<B->rmap->n; i++) { 4419 PetscInt row; 4420 row = i; 4421 if (rowindices) row = rowindices[i]; 4422 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4423 PetscInt col; 4424 col = Baij->j[count]; 4425 if (colindices) col = colindices[col]; 4426 v = Baij->a[count]; 4427 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4428 ++count; 4429 } 4430 } 4431 /* FIXME: set C's nonzerostate correctly. */ 4432 /* Assembly for C is necessary. */ 4433 C->preallocated = PETSC_TRUE; 4434 C->assembled = PETSC_TRUE; 4435 C->was_assembled = PETSC_FALSE; 4436 PetscFunctionReturn(0); 4437 } 4438 4439 4440 /* 4441 Special version for direct calls from Fortran 4442 */ 4443 #include <petsc/private/fortranimpl.h> 4444 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4445 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4446 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4447 #define matsetvaluesseqaij_ matsetvaluesseqaij 4448 #endif 4449 4450 /* Change these macros so can be used in void function */ 4451 #undef CHKERRQ 4452 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4453 #undef SETERRQ2 4454 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4455 #undef SETERRQ3 4456 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4457 4458 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) 4459 { 4460 Mat A = *AA; 4461 PetscInt m = *mm, n = *nn; 4462 InsertMode is = *isis; 4463 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4464 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4465 PetscInt *imax,*ai,*ailen; 4466 PetscErrorCode ierr; 4467 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4468 MatScalar *ap,value,*aa; 4469 PetscBool ignorezeroentries = a->ignorezeroentries; 4470 PetscBool roworiented = a->roworiented; 4471 4472 PetscFunctionBegin; 4473 MatCheckPreallocated(A,1); 4474 imax = a->imax; 4475 ai = a->i; 4476 ailen = a->ilen; 4477 aj = a->j; 4478 aa = a->a; 4479 4480 for (k=0; k<m; k++) { /* loop over added rows */ 4481 row = im[k]; 4482 if (row < 0) continue; 4483 #if defined(PETSC_USE_DEBUG) 4484 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4485 #endif 4486 rp = aj + ai[row]; ap = aa + ai[row]; 4487 rmax = imax[row]; nrow = ailen[row]; 4488 low = 0; 4489 high = nrow; 4490 for (l=0; l<n; l++) { /* loop over added columns */ 4491 if (in[l] < 0) continue; 4492 #if defined(PETSC_USE_DEBUG) 4493 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4494 #endif 4495 col = in[l]; 4496 if (roworiented) value = v[l + k*n]; 4497 else value = v[k + l*m]; 4498 4499 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4500 4501 if (col <= lastcol) low = 0; 4502 else high = nrow; 4503 lastcol = col; 4504 while (high-low > 5) { 4505 t = (low+high)/2; 4506 if (rp[t] > col) high = t; 4507 else low = t; 4508 } 4509 for (i=low; i<high; i++) { 4510 if (rp[i] > col) break; 4511 if (rp[i] == col) { 4512 if (is == ADD_VALUES) ap[i] += value; 4513 else ap[i] = value; 4514 goto noinsert; 4515 } 4516 } 4517 if (value == 0.0 && ignorezeroentries) goto noinsert; 4518 if (nonew == 1) goto noinsert; 4519 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4520 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4521 N = nrow++ - 1; a->nz++; high++; 4522 /* shift up all the later entries in this row */ 4523 for (ii=N; ii>=i; ii--) { 4524 rp[ii+1] = rp[ii]; 4525 ap[ii+1] = ap[ii]; 4526 } 4527 rp[i] = col; 4528 ap[i] = value; 4529 A->nonzerostate++; 4530 noinsert:; 4531 low = i + 1; 4532 } 4533 ailen[row] = nrow; 4534 } 4535 PetscFunctionReturnVoid(); 4536 } 4537 4538