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