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