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