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