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 case MAT_STRUCTURE_ONLY: 1109 /* These options are handled directly by MatSetOption() */ 1110 break; 1111 case MAT_NEW_DIAGONALS: 1112 case MAT_IGNORE_OFF_PROC_ENTRIES: 1113 case MAT_USE_HASH_TABLE: 1114 ierr = PetscInfo1(A,"Option %s ignored\n",MatOptions[op]);CHKERRQ(ierr); 1115 break; 1116 case MAT_USE_INODES: 1117 /* Not an error because MatSetOption_SeqAIJ_Inode handles this one */ 1118 break; 1119 case MAT_SUBMAT_SINGLEIS: 1120 A->submat_singleis = flg; 1121 break; 1122 default: 1123 SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_SUP,"unknown option %d",op); 1124 } 1125 ierr = MatSetOption_SeqAIJ_Inode(A,op,flg);CHKERRQ(ierr); 1126 PetscFunctionReturn(0); 1127 } 1128 1129 PetscErrorCode MatGetDiagonal_SeqAIJ(Mat A,Vec v) 1130 { 1131 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1132 PetscErrorCode ierr; 1133 PetscInt i,j,n,*ai=a->i,*aj=a->j,nz; 1134 PetscScalar *aa=a->a,*x,zero=0.0; 1135 1136 PetscFunctionBegin; 1137 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 1138 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 1139 1140 if (A->factortype == MAT_FACTOR_ILU || A->factortype == MAT_FACTOR_LU) { 1141 PetscInt *diag=a->diag; 1142 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1143 for (i=0; i<n; i++) x[i] = 1.0/aa[diag[i]]; 1144 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1145 PetscFunctionReturn(0); 1146 } 1147 1148 ierr = VecSet(v,zero);CHKERRQ(ierr); 1149 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 1150 for (i=0; i<n; i++) { 1151 nz = ai[i+1] - ai[i]; 1152 if (!nz) x[i] = 0.0; 1153 for (j=ai[i]; j<ai[i+1]; j++) { 1154 if (aj[j] == i) { 1155 x[i] = aa[j]; 1156 break; 1157 } 1158 } 1159 } 1160 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 1161 PetscFunctionReturn(0); 1162 } 1163 1164 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1165 PetscErrorCode MatMultTransposeAdd_SeqAIJ(Mat A,Vec xx,Vec zz,Vec yy) 1166 { 1167 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1168 PetscScalar *y; 1169 const PetscScalar *x; 1170 PetscErrorCode ierr; 1171 PetscInt m = A->rmap->n; 1172 #if !defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1173 const MatScalar *v; 1174 PetscScalar alpha; 1175 PetscInt n,i,j; 1176 const PetscInt *idx,*ii,*ridx=NULL; 1177 Mat_CompressedRow cprow = a->compressedrow; 1178 PetscBool usecprow = cprow.use; 1179 #endif 1180 1181 PetscFunctionBegin; 1182 if (zz != yy) {ierr = VecCopy(zz,yy);CHKERRQ(ierr);} 1183 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1184 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1185 1186 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTTRANSPOSEAIJ) 1187 fortranmulttransposeaddaij_(&m,x,a->i,a->j,a->a,y); 1188 #else 1189 if (usecprow) { 1190 m = cprow.nrows; 1191 ii = cprow.i; 1192 ridx = cprow.rindex; 1193 } else { 1194 ii = a->i; 1195 } 1196 for (i=0; i<m; i++) { 1197 idx = a->j + ii[i]; 1198 v = a->a + ii[i]; 1199 n = ii[i+1] - ii[i]; 1200 if (usecprow) { 1201 alpha = x[ridx[i]]; 1202 } else { 1203 alpha = x[i]; 1204 } 1205 for (j=0; j<n; j++) y[idx[j]] += alpha*v[j]; 1206 } 1207 #endif 1208 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1209 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1210 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1211 PetscFunctionReturn(0); 1212 } 1213 1214 PetscErrorCode MatMultTranspose_SeqAIJ(Mat A,Vec xx,Vec yy) 1215 { 1216 PetscErrorCode ierr; 1217 1218 PetscFunctionBegin; 1219 ierr = VecSet(yy,0.0);CHKERRQ(ierr); 1220 ierr = MatMultTransposeAdd_SeqAIJ(A,xx,yy,yy);CHKERRQ(ierr); 1221 PetscFunctionReturn(0); 1222 } 1223 1224 #include <../src/mat/impls/aij/seq/ftn-kernels/fmult.h> 1225 1226 PetscErrorCode MatMult_SeqAIJ(Mat A,Vec xx,Vec yy) 1227 { 1228 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1229 PetscScalar *y; 1230 const PetscScalar *x; 1231 const MatScalar *aa; 1232 PetscErrorCode ierr; 1233 PetscInt m=A->rmap->n; 1234 const PetscInt *aj,*ii,*ridx=NULL; 1235 PetscInt n,i; 1236 PetscScalar sum; 1237 PetscBool usecprow=a->compressedrow.use; 1238 1239 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1240 #pragma disjoint(*x,*y,*aa) 1241 #endif 1242 1243 PetscFunctionBegin; 1244 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1245 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1246 ii = a->i; 1247 if (usecprow) { /* use compressed row format */ 1248 ierr = PetscMemzero(y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1249 m = a->compressedrow.nrows; 1250 ii = a->compressedrow.i; 1251 ridx = a->compressedrow.rindex; 1252 for (i=0; i<m; i++) { 1253 n = ii[i+1] - ii[i]; 1254 aj = a->j + ii[i]; 1255 aa = a->a + ii[i]; 1256 sum = 0.0; 1257 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1258 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1259 y[*ridx++] = sum; 1260 } 1261 } else { /* do not use compressed row format */ 1262 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTAIJ) 1263 aj = a->j; 1264 aa = a->a; 1265 fortranmultaij_(&m,x,ii,aj,aa,y); 1266 #else 1267 for (i=0; i<m; i++) { 1268 n = ii[i+1] - ii[i]; 1269 aj = a->j + ii[i]; 1270 aa = a->a + ii[i]; 1271 sum = 0.0; 1272 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1273 y[i] = sum; 1274 } 1275 #endif 1276 } 1277 ierr = PetscLogFlops(2.0*a->nz - a->nonzerorowcnt);CHKERRQ(ierr); 1278 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1279 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1280 PetscFunctionReturn(0); 1281 } 1282 1283 PetscErrorCode MatMultMax_SeqAIJ(Mat A,Vec xx,Vec yy) 1284 { 1285 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1286 PetscScalar *y; 1287 const PetscScalar *x; 1288 const MatScalar *aa; 1289 PetscErrorCode ierr; 1290 PetscInt m=A->rmap->n; 1291 const PetscInt *aj,*ii,*ridx=NULL; 1292 PetscInt n,i,nonzerorow=0; 1293 PetscScalar sum; 1294 PetscBool usecprow=a->compressedrow.use; 1295 1296 #if defined(PETSC_HAVE_PRAGMA_DISJOINT) 1297 #pragma disjoint(*x,*y,*aa) 1298 #endif 1299 1300 PetscFunctionBegin; 1301 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1302 ierr = VecGetArray(yy,&y);CHKERRQ(ierr); 1303 if (usecprow) { /* use compressed row format */ 1304 m = a->compressedrow.nrows; 1305 ii = a->compressedrow.i; 1306 ridx = a->compressedrow.rindex; 1307 for (i=0; i<m; i++) { 1308 n = ii[i+1] - ii[i]; 1309 aj = a->j + ii[i]; 1310 aa = a->a + ii[i]; 1311 sum = 0.0; 1312 nonzerorow += (n>0); 1313 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1314 /* for (j=0; j<n; j++) sum += (*aa++)*x[*aj++]; */ 1315 y[*ridx++] = sum; 1316 } 1317 } else { /* do not use compressed row format */ 1318 ii = a->i; 1319 for (i=0; i<m; i++) { 1320 n = ii[i+1] - ii[i]; 1321 aj = a->j + ii[i]; 1322 aa = a->a + ii[i]; 1323 sum = 0.0; 1324 nonzerorow += (n>0); 1325 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1326 y[i] = sum; 1327 } 1328 } 1329 ierr = PetscLogFlops(2.0*a->nz - nonzerorow);CHKERRQ(ierr); 1330 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1331 ierr = VecRestoreArray(yy,&y);CHKERRQ(ierr); 1332 PetscFunctionReturn(0); 1333 } 1334 1335 PetscErrorCode MatMultAddMax_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1336 { 1337 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1338 PetscScalar *y,*z; 1339 const PetscScalar *x; 1340 const MatScalar *aa; 1341 PetscErrorCode ierr; 1342 PetscInt m = A->rmap->n,*aj,*ii; 1343 PetscInt n,i,*ridx=NULL; 1344 PetscScalar sum; 1345 PetscBool usecprow=a->compressedrow.use; 1346 1347 PetscFunctionBegin; 1348 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1349 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1350 if (usecprow) { /* use compressed row format */ 1351 if (zz != yy) { 1352 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1353 } 1354 m = a->compressedrow.nrows; 1355 ii = a->compressedrow.i; 1356 ridx = a->compressedrow.rindex; 1357 for (i=0; i<m; i++) { 1358 n = ii[i+1] - ii[i]; 1359 aj = a->j + ii[i]; 1360 aa = a->a + ii[i]; 1361 sum = y[*ridx]; 1362 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1363 z[*ridx++] = sum; 1364 } 1365 } else { /* do not use compressed row format */ 1366 ii = a->i; 1367 for (i=0; i<m; i++) { 1368 n = ii[i+1] - ii[i]; 1369 aj = a->j + ii[i]; 1370 aa = a->a + ii[i]; 1371 sum = y[i]; 1372 PetscSparseDenseMaxDot(sum,x,aa,aj,n); 1373 z[i] = sum; 1374 } 1375 } 1376 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1377 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1378 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1379 PetscFunctionReturn(0); 1380 } 1381 1382 #include <../src/mat/impls/aij/seq/ftn-kernels/fmultadd.h> 1383 PetscErrorCode MatMultAdd_SeqAIJ(Mat A,Vec xx,Vec yy,Vec zz) 1384 { 1385 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1386 PetscScalar *y,*z; 1387 const PetscScalar *x; 1388 const MatScalar *aa; 1389 PetscErrorCode ierr; 1390 const PetscInt *aj,*ii,*ridx=NULL; 1391 PetscInt m = A->rmap->n,n,i; 1392 PetscScalar sum; 1393 PetscBool usecprow=a->compressedrow.use; 1394 1395 PetscFunctionBegin; 1396 ierr = VecGetArrayRead(xx,&x);CHKERRQ(ierr); 1397 ierr = VecGetArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1398 if (usecprow) { /* use compressed row format */ 1399 if (zz != yy) { 1400 ierr = PetscMemcpy(z,y,m*sizeof(PetscScalar));CHKERRQ(ierr); 1401 } 1402 m = a->compressedrow.nrows; 1403 ii = a->compressedrow.i; 1404 ridx = a->compressedrow.rindex; 1405 for (i=0; i<m; i++) { 1406 n = ii[i+1] - ii[i]; 1407 aj = a->j + ii[i]; 1408 aa = a->a + ii[i]; 1409 sum = y[*ridx]; 1410 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1411 z[*ridx++] = sum; 1412 } 1413 } else { /* do not use compressed row format */ 1414 ii = a->i; 1415 #if defined(PETSC_USE_FORTRAN_KERNEL_MULTADDAIJ) 1416 aj = a->j; 1417 aa = a->a; 1418 fortranmultaddaij_(&m,x,ii,aj,aa,y,z); 1419 #else 1420 for (i=0; i<m; i++) { 1421 n = ii[i+1] - ii[i]; 1422 aj = a->j + ii[i]; 1423 aa = a->a + ii[i]; 1424 sum = y[i]; 1425 PetscSparseDensePlusDot(sum,x,aa,aj,n); 1426 z[i] = sum; 1427 } 1428 #endif 1429 } 1430 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1431 ierr = VecRestoreArrayRead(xx,&x);CHKERRQ(ierr); 1432 ierr = VecRestoreArrayPair(yy,zz,&y,&z);CHKERRQ(ierr); 1433 #if defined(PETSC_HAVE_CUSP) 1434 /* 1435 ierr = VecView(xx,0);CHKERRQ(ierr); 1436 ierr = VecView(zz,0);CHKERRQ(ierr); 1437 ierr = MatView(A,0);CHKERRQ(ierr); 1438 */ 1439 #endif 1440 PetscFunctionReturn(0); 1441 } 1442 1443 /* 1444 Adds diagonal pointers to sparse matrix structure. 1445 */ 1446 PetscErrorCode MatMarkDiagonal_SeqAIJ(Mat A) 1447 { 1448 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1449 PetscErrorCode ierr; 1450 PetscInt i,j,m = A->rmap->n; 1451 1452 PetscFunctionBegin; 1453 if (!a->diag) { 1454 ierr = PetscMalloc1(m,&a->diag);CHKERRQ(ierr); 1455 ierr = PetscLogObjectMemory((PetscObject)A, m*sizeof(PetscInt));CHKERRQ(ierr); 1456 } 1457 for (i=0; i<A->rmap->n; i++) { 1458 a->diag[i] = a->i[i+1]; 1459 for (j=a->i[i]; j<a->i[i+1]; j++) { 1460 if (a->j[j] == i) { 1461 a->diag[i] = j; 1462 break; 1463 } 1464 } 1465 } 1466 PetscFunctionReturn(0); 1467 } 1468 1469 /* 1470 Checks for missing diagonals 1471 */ 1472 PetscErrorCode MatMissingDiagonal_SeqAIJ(Mat A,PetscBool *missing,PetscInt *d) 1473 { 1474 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1475 PetscInt *diag,*ii = a->i,i; 1476 1477 PetscFunctionBegin; 1478 *missing = PETSC_FALSE; 1479 if (A->rmap->n > 0 && !ii) { 1480 *missing = PETSC_TRUE; 1481 if (d) *d = 0; 1482 PetscInfo(A,"Matrix has no entries therefore is missing diagonal\n"); 1483 } else { 1484 diag = a->diag; 1485 for (i=0; i<A->rmap->n; i++) { 1486 if (diag[i] >= ii[i+1]) { 1487 *missing = PETSC_TRUE; 1488 if (d) *d = i; 1489 PetscInfo1(A,"Matrix is missing diagonal number %D\n",i); 1490 break; 1491 } 1492 } 1493 } 1494 PetscFunctionReturn(0); 1495 } 1496 1497 /* 1498 Negative shift indicates do not generate an error if there is a zero diagonal, just invert it anyways 1499 */ 1500 PetscErrorCode MatInvertDiagonal_SeqAIJ(Mat A,PetscScalar omega,PetscScalar fshift) 1501 { 1502 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 1503 PetscErrorCode ierr; 1504 PetscInt i,*diag,m = A->rmap->n; 1505 MatScalar *v = a->a; 1506 PetscScalar *idiag,*mdiag; 1507 1508 PetscFunctionBegin; 1509 if (a->idiagvalid) PetscFunctionReturn(0); 1510 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 1511 diag = a->diag; 1512 if (!a->idiag) { 1513 ierr = PetscMalloc3(m,&a->idiag,m,&a->mdiag,m,&a->ssor_work);CHKERRQ(ierr); 1514 ierr = PetscLogObjectMemory((PetscObject)A, 3*m*sizeof(PetscScalar));CHKERRQ(ierr); 1515 v = a->a; 1516 } 1517 mdiag = a->mdiag; 1518 idiag = a->idiag; 1519 1520 if (omega == 1.0 && PetscRealPart(fshift) <= 0.0) { 1521 for (i=0; i<m; i++) { 1522 mdiag[i] = v[diag[i]]; 1523 if (!PetscAbsScalar(mdiag[i])) { /* zero diagonal */ 1524 if (PetscRealPart(fshift)) { 1525 ierr = PetscInfo1(A,"Zero diagonal on row %D\n",i);CHKERRQ(ierr); 1526 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 1527 A->factorerror_zeropivot_value = 0.0; 1528 A->factorerror_zeropivot_row = i; 1529 } SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Zero diagonal on row %D",i); 1530 } 1531 idiag[i] = 1.0/v[diag[i]]; 1532 } 1533 ierr = PetscLogFlops(m);CHKERRQ(ierr); 1534 } else { 1535 for (i=0; i<m; i++) { 1536 mdiag[i] = v[diag[i]]; 1537 idiag[i] = omega/(fshift + v[diag[i]]); 1538 } 1539 ierr = PetscLogFlops(2.0*m);CHKERRQ(ierr); 1540 } 1541 a->idiagvalid = PETSC_TRUE; 1542 PetscFunctionReturn(0); 1543 } 1544 1545 #include <../src/mat/impls/aij/seq/ftn-kernels/frelax.h> 1546 PetscErrorCode MatSOR_SeqAIJ(Mat A,Vec bb,PetscReal omega,MatSORType flag,PetscReal fshift,PetscInt its,PetscInt lits,Vec xx) 1547 { 1548 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1549 PetscScalar *x,d,sum,*t,scale; 1550 const MatScalar *v,*idiag=0,*mdiag; 1551 const PetscScalar *b, *bs,*xb, *ts; 1552 PetscErrorCode ierr; 1553 PetscInt n,m = A->rmap->n,i; 1554 const PetscInt *idx,*diag; 1555 1556 PetscFunctionBegin; 1557 its = its*lits; 1558 1559 if (fshift != a->fshift || omega != a->omega) a->idiagvalid = PETSC_FALSE; /* must recompute idiag[] */ 1560 if (!a->idiagvalid) {ierr = MatInvertDiagonal_SeqAIJ(A,omega,fshift);CHKERRQ(ierr);} 1561 a->fshift = fshift; 1562 a->omega = omega; 1563 1564 diag = a->diag; 1565 t = a->ssor_work; 1566 idiag = a->idiag; 1567 mdiag = a->mdiag; 1568 1569 ierr = VecGetArray(xx,&x);CHKERRQ(ierr); 1570 ierr = VecGetArrayRead(bb,&b);CHKERRQ(ierr); 1571 /* We count flops by assuming the upper triangular and lower triangular parts have the same number of nonzeros */ 1572 if (flag == SOR_APPLY_UPPER) { 1573 /* apply (U + D/omega) to the vector */ 1574 bs = b; 1575 for (i=0; i<m; i++) { 1576 d = fshift + mdiag[i]; 1577 n = a->i[i+1] - diag[i] - 1; 1578 idx = a->j + diag[i] + 1; 1579 v = a->a + diag[i] + 1; 1580 sum = b[i]*d/omega; 1581 PetscSparseDensePlusDot(sum,bs,v,idx,n); 1582 x[i] = sum; 1583 } 1584 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1585 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1586 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1587 PetscFunctionReturn(0); 1588 } 1589 1590 if (flag == SOR_APPLY_LOWER) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"SOR_APPLY_LOWER is not implemented"); 1591 else if (flag & SOR_EISENSTAT) { 1592 /* Let A = L + U + D; where L is lower trianglar, 1593 U is upper triangular, E = D/omega; This routine applies 1594 1595 (L + E)^{-1} A (U + E)^{-1} 1596 1597 to a vector efficiently using Eisenstat's trick. 1598 */ 1599 scale = (2.0/omega) - 1.0; 1600 1601 /* x = (E + U)^{-1} b */ 1602 for (i=m-1; i>=0; i--) { 1603 n = a->i[i+1] - diag[i] - 1; 1604 idx = a->j + diag[i] + 1; 1605 v = a->a + diag[i] + 1; 1606 sum = b[i]; 1607 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1608 x[i] = sum*idiag[i]; 1609 } 1610 1611 /* t = b - (2*E - D)x */ 1612 v = a->a; 1613 for (i=0; i<m; i++) t[i] = b[i] - scale*(v[*diag++])*x[i]; 1614 1615 /* t = (E + L)^{-1}t */ 1616 ts = t; 1617 diag = a->diag; 1618 for (i=0; i<m; i++) { 1619 n = diag[i] - a->i[i]; 1620 idx = a->j + a->i[i]; 1621 v = a->a + a->i[i]; 1622 sum = t[i]; 1623 PetscSparseDenseMinusDot(sum,ts,v,idx,n); 1624 t[i] = sum*idiag[i]; 1625 /* x = x + t */ 1626 x[i] += t[i]; 1627 } 1628 1629 ierr = PetscLogFlops(6.0*m-1 + 2.0*a->nz);CHKERRQ(ierr); 1630 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1631 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1632 PetscFunctionReturn(0); 1633 } 1634 if (flag & SOR_ZERO_INITIAL_GUESS) { 1635 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1636 for (i=0; i<m; i++) { 1637 n = diag[i] - a->i[i]; 1638 idx = a->j + a->i[i]; 1639 v = a->a + a->i[i]; 1640 sum = b[i]; 1641 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1642 t[i] = sum; 1643 x[i] = sum*idiag[i]; 1644 } 1645 xb = t; 1646 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 1647 } else xb = b; 1648 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1649 for (i=m-1; i>=0; i--) { 1650 n = a->i[i+1] - diag[i] - 1; 1651 idx = a->j + diag[i] + 1; 1652 v = a->a + diag[i] + 1; 1653 sum = xb[i]; 1654 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1655 if (xb == b) { 1656 x[i] = sum*idiag[i]; 1657 } else { 1658 x[i] = (1-omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1659 } 1660 } 1661 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1662 } 1663 its--; 1664 } 1665 while (its--) { 1666 if (flag & SOR_FORWARD_SWEEP || flag & SOR_LOCAL_FORWARD_SWEEP) { 1667 for (i=0; i<m; i++) { 1668 /* lower */ 1669 n = diag[i] - a->i[i]; 1670 idx = a->j + a->i[i]; 1671 v = a->a + a->i[i]; 1672 sum = b[i]; 1673 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1674 t[i] = sum; /* save application of the lower-triangular part */ 1675 /* upper */ 1676 n = a->i[i+1] - diag[i] - 1; 1677 idx = a->j + diag[i] + 1; 1678 v = a->a + diag[i] + 1; 1679 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1680 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1681 } 1682 xb = t; 1683 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1684 } else xb = b; 1685 if (flag & SOR_BACKWARD_SWEEP || flag & SOR_LOCAL_BACKWARD_SWEEP) { 1686 for (i=m-1; i>=0; i--) { 1687 sum = xb[i]; 1688 if (xb == b) { 1689 /* whole matrix (no checkpointing available) */ 1690 n = a->i[i+1] - a->i[i]; 1691 idx = a->j + a->i[i]; 1692 v = a->a + a->i[i]; 1693 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1694 x[i] = (1. - omega)*x[i] + (sum + mdiag[i]*x[i])*idiag[i]; 1695 } else { /* lower-triangular part has been saved, so only apply upper-triangular */ 1696 n = a->i[i+1] - diag[i] - 1; 1697 idx = a->j + diag[i] + 1; 1698 v = a->a + diag[i] + 1; 1699 PetscSparseDenseMinusDot(sum,x,v,idx,n); 1700 x[i] = (1. - omega)*x[i] + sum*idiag[i]; /* omega in idiag */ 1701 } 1702 } 1703 if (xb == b) { 1704 ierr = PetscLogFlops(2.0*a->nz);CHKERRQ(ierr); 1705 } else { 1706 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); /* assumes 1/2 in upper */ 1707 } 1708 } 1709 } 1710 ierr = VecRestoreArray(xx,&x);CHKERRQ(ierr); 1711 ierr = VecRestoreArrayRead(bb,&b);CHKERRQ(ierr); 1712 PetscFunctionReturn(0); 1713 } 1714 1715 1716 PetscErrorCode MatGetInfo_SeqAIJ(Mat A,MatInfoType flag,MatInfo *info) 1717 { 1718 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1719 1720 PetscFunctionBegin; 1721 info->block_size = 1.0; 1722 info->nz_allocated = (double)a->maxnz; 1723 info->nz_used = (double)a->nz; 1724 info->nz_unneeded = (double)(a->maxnz - a->nz); 1725 info->assemblies = (double)A->num_ass; 1726 info->mallocs = (double)A->info.mallocs; 1727 info->memory = ((PetscObject)A)->mem; 1728 if (A->factortype) { 1729 info->fill_ratio_given = A->info.fill_ratio_given; 1730 info->fill_ratio_needed = A->info.fill_ratio_needed; 1731 info->factor_mallocs = A->info.factor_mallocs; 1732 } else { 1733 info->fill_ratio_given = 0; 1734 info->fill_ratio_needed = 0; 1735 info->factor_mallocs = 0; 1736 } 1737 PetscFunctionReturn(0); 1738 } 1739 1740 PetscErrorCode MatZeroRows_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1741 { 1742 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1743 PetscInt i,m = A->rmap->n - 1; 1744 PetscErrorCode ierr; 1745 const PetscScalar *xx; 1746 PetscScalar *bb; 1747 PetscInt d = 0; 1748 1749 PetscFunctionBegin; 1750 if (x && b) { 1751 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1752 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1753 for (i=0; i<N; i++) { 1754 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1755 bb[rows[i]] = diag*xx[rows[i]]; 1756 } 1757 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1758 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1759 } 1760 1761 if (a->keepnonzeropattern) { 1762 for (i=0; i<N; i++) { 1763 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1764 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1765 } 1766 if (diag != 0.0) { 1767 for (i=0; i<N; i++) { 1768 d = rows[i]; 1769 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); 1770 } 1771 for (i=0; i<N; i++) { 1772 a->a[a->diag[rows[i]]] = diag; 1773 } 1774 } 1775 } else { 1776 if (diag != 0.0) { 1777 for (i=0; i<N; i++) { 1778 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1779 if (a->ilen[rows[i]] > 0) { 1780 a->ilen[rows[i]] = 1; 1781 a->a[a->i[rows[i]]] = diag; 1782 a->j[a->i[rows[i]]] = rows[i]; 1783 } else { /* in case row was completely empty */ 1784 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1785 } 1786 } 1787 } else { 1788 for (i=0; i<N; i++) { 1789 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1790 a->ilen[rows[i]] = 0; 1791 } 1792 } 1793 A->nonzerostate++; 1794 } 1795 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1796 PetscFunctionReturn(0); 1797 } 1798 1799 PetscErrorCode MatZeroRowsColumns_SeqAIJ(Mat A,PetscInt N,const PetscInt rows[],PetscScalar diag,Vec x,Vec b) 1800 { 1801 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1802 PetscInt i,j,m = A->rmap->n - 1,d = 0; 1803 PetscErrorCode ierr; 1804 PetscBool missing,*zeroed,vecs = PETSC_FALSE; 1805 const PetscScalar *xx; 1806 PetscScalar *bb; 1807 1808 PetscFunctionBegin; 1809 if (x && b) { 1810 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr); 1811 ierr = VecGetArray(b,&bb);CHKERRQ(ierr); 1812 vecs = PETSC_TRUE; 1813 } 1814 ierr = PetscCalloc1(A->rmap->n,&zeroed);CHKERRQ(ierr); 1815 for (i=0; i<N; i++) { 1816 if (rows[i] < 0 || rows[i] > m) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"row %D out of range", rows[i]); 1817 ierr = PetscMemzero(&a->a[a->i[rows[i]]],a->ilen[rows[i]]*sizeof(PetscScalar));CHKERRQ(ierr); 1818 1819 zeroed[rows[i]] = PETSC_TRUE; 1820 } 1821 for (i=0; i<A->rmap->n; i++) { 1822 if (!zeroed[i]) { 1823 for (j=a->i[i]; j<a->i[i+1]; j++) { 1824 if (zeroed[a->j[j]]) { 1825 if (vecs) bb[i] -= a->a[j]*xx[a->j[j]]; 1826 a->a[j] = 0.0; 1827 } 1828 } 1829 } else if (vecs) bb[i] = diag*xx[i]; 1830 } 1831 if (x && b) { 1832 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr); 1833 ierr = VecRestoreArray(b,&bb);CHKERRQ(ierr); 1834 } 1835 ierr = PetscFree(zeroed);CHKERRQ(ierr); 1836 if (diag != 0.0) { 1837 ierr = MatMissingDiagonal_SeqAIJ(A,&missing,&d);CHKERRQ(ierr); 1838 if (missing) { 1839 if (a->nonew) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Matrix is missing diagonal entry in row %D",d); 1840 else { 1841 for (i=0; i<N; i++) { 1842 ierr = MatSetValues_SeqAIJ(A,1,&rows[i],1,&rows[i],&diag,INSERT_VALUES);CHKERRQ(ierr); 1843 } 1844 } 1845 } else { 1846 for (i=0; i<N; i++) { 1847 a->a[a->diag[rows[i]]] = diag; 1848 } 1849 } 1850 } 1851 ierr = MatAssemblyEnd_SeqAIJ(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1852 PetscFunctionReturn(0); 1853 } 1854 1855 PetscErrorCode MatGetRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1856 { 1857 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1858 PetscInt *itmp; 1859 1860 PetscFunctionBegin; 1861 if (row < 0 || row >= A->rmap->n) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Row %D out of range",row); 1862 1863 *nz = a->i[row+1] - a->i[row]; 1864 if (v) *v = a->a + a->i[row]; 1865 if (idx) { 1866 itmp = a->j + a->i[row]; 1867 if (*nz) *idx = itmp; 1868 else *idx = 0; 1869 } 1870 PetscFunctionReturn(0); 1871 } 1872 1873 /* remove this function? */ 1874 PetscErrorCode MatRestoreRow_SeqAIJ(Mat A,PetscInt row,PetscInt *nz,PetscInt **idx,PetscScalar **v) 1875 { 1876 PetscFunctionBegin; 1877 PetscFunctionReturn(0); 1878 } 1879 1880 PetscErrorCode MatNorm_SeqAIJ(Mat A,NormType type,PetscReal *nrm) 1881 { 1882 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1883 MatScalar *v = a->a; 1884 PetscReal sum = 0.0; 1885 PetscErrorCode ierr; 1886 PetscInt i,j; 1887 1888 PetscFunctionBegin; 1889 if (type == NORM_FROBENIUS) { 1890 #if defined(PETSC_USE_REAL___FP16) 1891 PetscBLASInt one = 1,nz = a->nz; 1892 *nrm = BLASnrm2_(&nz,v,&one); 1893 #else 1894 for (i=0; i<a->nz; i++) { 1895 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1896 } 1897 *nrm = PetscSqrtReal(sum); 1898 #endif 1899 ierr = PetscLogFlops(2*a->nz);CHKERRQ(ierr); 1900 } else if (type == NORM_1) { 1901 PetscReal *tmp; 1902 PetscInt *jj = a->j; 1903 ierr = PetscCalloc1(A->cmap->n+1,&tmp);CHKERRQ(ierr); 1904 *nrm = 0.0; 1905 for (j=0; j<a->nz; j++) { 1906 tmp[*jj++] += PetscAbsScalar(*v); v++; 1907 } 1908 for (j=0; j<A->cmap->n; j++) { 1909 if (tmp[j] > *nrm) *nrm = tmp[j]; 1910 } 1911 ierr = PetscFree(tmp);CHKERRQ(ierr); 1912 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1913 } else if (type == NORM_INFINITY) { 1914 *nrm = 0.0; 1915 for (j=0; j<A->rmap->n; j++) { 1916 v = a->a + a->i[j]; 1917 sum = 0.0; 1918 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1919 sum += PetscAbsScalar(*v); v++; 1920 } 1921 if (sum > *nrm) *nrm = sum; 1922 } 1923 ierr = PetscLogFlops(PetscMax(a->nz-1,0));CHKERRQ(ierr); 1924 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1925 PetscFunctionReturn(0); 1926 } 1927 1928 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1929 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1930 { 1931 PetscErrorCode ierr; 1932 PetscInt i,j,anzj; 1933 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data,*b; 1934 PetscInt an=A->cmap->N,am=A->rmap->N; 1935 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1936 1937 PetscFunctionBegin; 1938 /* Allocate space for symbolic transpose info and work array */ 1939 ierr = PetscCalloc1(an+1,&ati);CHKERRQ(ierr); 1940 ierr = PetscMalloc1(ai[am],&atj);CHKERRQ(ierr); 1941 ierr = PetscMalloc1(an,&atfill);CHKERRQ(ierr); 1942 1943 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 1944 /* Note: offset by 1 for fast conversion into csr format. */ 1945 for (i=0;i<ai[am];i++) ati[aj[i]+1] += 1; 1946 /* Form ati for csr format of A^T. */ 1947 for (i=0;i<an;i++) ati[i+1] += ati[i]; 1948 1949 /* Copy ati into atfill so we have locations of the next free space in atj */ 1950 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 1951 1952 /* Walk through A row-wise and mark nonzero entries of A^T. */ 1953 for (i=0;i<am;i++) { 1954 anzj = ai[i+1] - ai[i]; 1955 for (j=0;j<anzj;j++) { 1956 atj[atfill[*aj]] = i; 1957 atfill[*aj++] += 1; 1958 } 1959 } 1960 1961 /* Clean up temporary space and complete requests. */ 1962 ierr = PetscFree(atfill);CHKERRQ(ierr); 1963 ierr = MatCreateSeqAIJWithArrays(PetscObjectComm((PetscObject)A),an,am,ati,atj,NULL,B);CHKERRQ(ierr); 1964 ierr = MatSetBlockSizes(*B,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 1965 1966 b = (Mat_SeqAIJ*)((*B)->data); 1967 b->free_a = PETSC_FALSE; 1968 b->free_ij = PETSC_TRUE; 1969 b->nonew = 0; 1970 PetscFunctionReturn(0); 1971 } 1972 1973 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1974 { 1975 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1976 Mat C; 1977 PetscErrorCode ierr; 1978 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1979 MatScalar *array = a->a; 1980 1981 PetscFunctionBegin; 1982 if (reuse == MAT_INPLACE_MATRIX && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1983 1984 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_INPLACE_MATRIX) { 1985 ierr = PetscCalloc1(1+A->cmap->n,&col);CHKERRQ(ierr); 1986 1987 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1988 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 1989 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1990 ierr = MatSetBlockSizes(C,PetscAbs(A->cmap->bs),PetscAbs(A->rmap->bs));CHKERRQ(ierr); 1991 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1992 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1993 ierr = PetscFree(col);CHKERRQ(ierr); 1994 } else { 1995 C = *B; 1996 } 1997 1998 for (i=0; i<m; i++) { 1999 len = ai[i+1]-ai[i]; 2000 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 2001 array += len; 2002 aj += len; 2003 } 2004 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2005 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2006 2007 if (reuse == MAT_INITIAL_MATRIX || reuse == MAT_REUSE_MATRIX) { 2008 *B = C; 2009 } else { 2010 ierr = MatHeaderMerge(A,&C);CHKERRQ(ierr); 2011 } 2012 PetscFunctionReturn(0); 2013 } 2014 2015 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2016 { 2017 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2018 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2019 MatScalar *va,*vb; 2020 PetscErrorCode ierr; 2021 PetscInt ma,na,mb,nb, i; 2022 2023 PetscFunctionBegin; 2024 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2025 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2026 if (ma!=nb || na!=mb) { 2027 *f = PETSC_FALSE; 2028 PetscFunctionReturn(0); 2029 } 2030 aii = aij->i; bii = bij->i; 2031 adx = aij->j; bdx = bij->j; 2032 va = aij->a; vb = bij->a; 2033 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2034 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2035 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2036 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2037 2038 *f = PETSC_TRUE; 2039 for (i=0; i<ma; i++) { 2040 while (aptr[i]<aii[i+1]) { 2041 PetscInt idc,idr; 2042 PetscScalar vc,vr; 2043 /* column/row index/value */ 2044 idc = adx[aptr[i]]; 2045 idr = bdx[bptr[idc]]; 2046 vc = va[aptr[i]]; 2047 vr = vb[bptr[idc]]; 2048 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 2049 *f = PETSC_FALSE; 2050 goto done; 2051 } else { 2052 aptr[i]++; 2053 if (B || i!=idc) bptr[idc]++; 2054 } 2055 } 2056 } 2057 done: 2058 ierr = PetscFree(aptr);CHKERRQ(ierr); 2059 ierr = PetscFree(bptr);CHKERRQ(ierr); 2060 PetscFunctionReturn(0); 2061 } 2062 2063 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 2064 { 2065 Mat_SeqAIJ *aij = (Mat_SeqAIJ*) A->data,*bij = (Mat_SeqAIJ*) B->data; 2066 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 2067 MatScalar *va,*vb; 2068 PetscErrorCode ierr; 2069 PetscInt ma,na,mb,nb, i; 2070 2071 PetscFunctionBegin; 2072 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2073 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2074 if (ma!=nb || na!=mb) { 2075 *f = PETSC_FALSE; 2076 PetscFunctionReturn(0); 2077 } 2078 aii = aij->i; bii = bij->i; 2079 adx = aij->j; bdx = bij->j; 2080 va = aij->a; vb = bij->a; 2081 ierr = PetscMalloc1(ma,&aptr);CHKERRQ(ierr); 2082 ierr = PetscMalloc1(mb,&bptr);CHKERRQ(ierr); 2083 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2084 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2085 2086 *f = PETSC_TRUE; 2087 for (i=0; i<ma; i++) { 2088 while (aptr[i]<aii[i+1]) { 2089 PetscInt idc,idr; 2090 PetscScalar vc,vr; 2091 /* column/row index/value */ 2092 idc = adx[aptr[i]]; 2093 idr = bdx[bptr[idc]]; 2094 vc = va[aptr[i]]; 2095 vr = vb[bptr[idc]]; 2096 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2097 *f = PETSC_FALSE; 2098 goto done; 2099 } else { 2100 aptr[i]++; 2101 if (B || i!=idc) bptr[idc]++; 2102 } 2103 } 2104 } 2105 done: 2106 ierr = PetscFree(aptr);CHKERRQ(ierr); 2107 ierr = PetscFree(bptr);CHKERRQ(ierr); 2108 PetscFunctionReturn(0); 2109 } 2110 2111 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2112 { 2113 PetscErrorCode ierr; 2114 2115 PetscFunctionBegin; 2116 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2117 PetscFunctionReturn(0); 2118 } 2119 2120 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2121 { 2122 PetscErrorCode ierr; 2123 2124 PetscFunctionBegin; 2125 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2126 PetscFunctionReturn(0); 2127 } 2128 2129 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2130 { 2131 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2132 PetscScalar *l,*r,x; 2133 MatScalar *v; 2134 PetscErrorCode ierr; 2135 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2136 2137 PetscFunctionBegin; 2138 if (ll) { 2139 /* The local size is used so that VecMPI can be passed to this routine 2140 by MatDiagonalScale_MPIAIJ */ 2141 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2142 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2143 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2144 v = a->a; 2145 for (i=0; i<m; i++) { 2146 x = l[i]; 2147 M = a->i[i+1] - a->i[i]; 2148 for (j=0; j<M; j++) (*v++) *= x; 2149 } 2150 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2151 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2152 } 2153 if (rr) { 2154 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2155 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2156 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2157 v = a->a; jj = a->j; 2158 for (i=0; i<nz; i++) (*v++) *= r[*jj++]; 2159 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2160 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2161 } 2162 ierr = MatSeqAIJInvalidateDiagonal(A);CHKERRQ(ierr); 2163 PetscFunctionReturn(0); 2164 } 2165 2166 PetscErrorCode MatCreateSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2167 { 2168 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2169 PetscErrorCode ierr; 2170 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2171 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2172 const PetscInt *irow,*icol; 2173 PetscInt nrows,ncols; 2174 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2175 MatScalar *a_new,*mat_a; 2176 Mat C; 2177 PetscBool stride; 2178 2179 PetscFunctionBegin; 2180 2181 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2182 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2183 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2184 2185 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2186 if (stride) { 2187 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2188 } else { 2189 first = 0; 2190 step = 0; 2191 } 2192 if (stride && step == 1) { 2193 /* special case of contiguous rows */ 2194 ierr = PetscMalloc2(nrows,&lens,nrows,&starts);CHKERRQ(ierr); 2195 /* loop over new rows determining lens and starting points */ 2196 for (i=0; i<nrows; i++) { 2197 kstart = ai[irow[i]]; 2198 kend = kstart + ailen[irow[i]]; 2199 starts[i] = kstart; 2200 for (k=kstart; k<kend; k++) { 2201 if (aj[k] >= first) { 2202 starts[i] = k; 2203 break; 2204 } 2205 } 2206 sum = 0; 2207 while (k < kend) { 2208 if (aj[k++] >= first+ncols) break; 2209 sum++; 2210 } 2211 lens[i] = sum; 2212 } 2213 /* create submatrix */ 2214 if (scall == MAT_REUSE_MATRIX) { 2215 PetscInt n_cols,n_rows; 2216 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2217 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2218 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2219 C = *B; 2220 } else { 2221 PetscInt rbs,cbs; 2222 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2223 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2224 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2225 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2226 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2227 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2228 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2229 } 2230 c = (Mat_SeqAIJ*)C->data; 2231 2232 /* loop over rows inserting into submatrix */ 2233 a_new = c->a; 2234 j_new = c->j; 2235 i_new = c->i; 2236 2237 for (i=0; i<nrows; i++) { 2238 ii = starts[i]; 2239 lensi = lens[i]; 2240 for (k=0; k<lensi; k++) { 2241 *j_new++ = aj[ii+k] - first; 2242 } 2243 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2244 a_new += lensi; 2245 i_new[i+1] = i_new[i] + lensi; 2246 c->ilen[i] = lensi; 2247 } 2248 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2249 } else { 2250 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2251 ierr = PetscCalloc1(oldcols,&smap);CHKERRQ(ierr); 2252 ierr = PetscMalloc1(1+nrows,&lens);CHKERRQ(ierr); 2253 for (i=0; i<ncols; i++) { 2254 #if defined(PETSC_USE_DEBUG) 2255 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); 2256 #endif 2257 smap[icol[i]] = i+1; 2258 } 2259 2260 /* determine lens of each row */ 2261 for (i=0; i<nrows; i++) { 2262 kstart = ai[irow[i]]; 2263 kend = kstart + a->ilen[irow[i]]; 2264 lens[i] = 0; 2265 for (k=kstart; k<kend; k++) { 2266 if (smap[aj[k]]) { 2267 lens[i]++; 2268 } 2269 } 2270 } 2271 /* Create and fill new matrix */ 2272 if (scall == MAT_REUSE_MATRIX) { 2273 PetscBool equal; 2274 2275 c = (Mat_SeqAIJ*)((*B)->data); 2276 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2277 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2278 if (!equal) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2279 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2280 C = *B; 2281 } else { 2282 PetscInt rbs,cbs; 2283 ierr = MatCreate(PetscObjectComm((PetscObject)A),&C);CHKERRQ(ierr); 2284 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2285 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2286 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2287 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2288 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2289 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2290 } 2291 c = (Mat_SeqAIJ*)(C->data); 2292 for (i=0; i<nrows; i++) { 2293 row = irow[i]; 2294 kstart = ai[row]; 2295 kend = kstart + a->ilen[row]; 2296 mat_i = c->i[i]; 2297 mat_j = c->j + mat_i; 2298 mat_a = c->a + mat_i; 2299 mat_ilen = c->ilen + i; 2300 for (k=kstart; k<kend; k++) { 2301 if ((tcol=smap[a->j[k]])) { 2302 *mat_j++ = tcol - 1; 2303 *mat_a++ = a->a[k]; 2304 (*mat_ilen)++; 2305 2306 } 2307 } 2308 } 2309 /* Free work space */ 2310 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2311 ierr = PetscFree(smap);CHKERRQ(ierr); 2312 ierr = PetscFree(lens);CHKERRQ(ierr); 2313 /* sort */ 2314 for (i = 0; i < nrows; i++) { 2315 PetscInt ilen; 2316 2317 mat_i = c->i[i]; 2318 mat_j = c->j + mat_i; 2319 mat_a = c->a + mat_i; 2320 ilen = c->ilen[i]; 2321 ierr = PetscSortIntWithScalarArray(ilen,mat_j,mat_a);CHKERRQ(ierr); 2322 } 2323 } 2324 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2325 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2326 2327 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2328 *B = C; 2329 PetscFunctionReturn(0); 2330 } 2331 2332 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat *subMat) 2333 { 2334 PetscErrorCode ierr; 2335 Mat B; 2336 2337 PetscFunctionBegin; 2338 if (scall == MAT_INITIAL_MATRIX) { 2339 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2340 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2341 ierr = MatSetBlockSizesFromMats(B,mat,mat);CHKERRQ(ierr); 2342 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2343 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2344 *subMat = B; 2345 } else { 2346 ierr = MatCopy_SeqAIJ(mat,*subMat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 2347 } 2348 PetscFunctionReturn(0); 2349 } 2350 2351 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2352 { 2353 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2354 PetscErrorCode ierr; 2355 Mat outA; 2356 PetscBool row_identity,col_identity; 2357 2358 PetscFunctionBegin; 2359 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2360 2361 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2362 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2363 2364 outA = inA; 2365 outA->factortype = MAT_FACTOR_LU; 2366 ierr = PetscFree(inA->solvertype);CHKERRQ(ierr); 2367 ierr = PetscStrallocpy(MATSOLVERPETSC,&inA->solvertype);CHKERRQ(ierr); 2368 2369 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2370 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2371 2372 a->row = row; 2373 2374 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2375 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2376 2377 a->col = col; 2378 2379 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2380 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2381 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2382 ierr = PetscLogObjectParent((PetscObject)inA,(PetscObject)a->icol);CHKERRQ(ierr); 2383 2384 if (!a->solve_work) { /* this matrix may have been factored before */ 2385 ierr = PetscMalloc1(inA->rmap->n+1,&a->solve_work);CHKERRQ(ierr); 2386 ierr = PetscLogObjectMemory((PetscObject)inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2387 } 2388 2389 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2390 if (row_identity && col_identity) { 2391 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2392 } else { 2393 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2394 } 2395 PetscFunctionReturn(0); 2396 } 2397 2398 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2399 { 2400 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2401 PetscScalar oalpha = alpha; 2402 PetscErrorCode ierr; 2403 PetscBLASInt one = 1,bnz; 2404 2405 PetscFunctionBegin; 2406 ierr = PetscBLASIntCast(a->nz,&bnz);CHKERRQ(ierr); 2407 PetscStackCallBLAS("BLASscal",BLASscal_(&bnz,&oalpha,a->a,&one)); 2408 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2409 ierr = MatSeqAIJInvalidateDiagonal(inA);CHKERRQ(ierr); 2410 PetscFunctionReturn(0); 2411 } 2412 2413 PetscErrorCode MatDestroySubMatrices_Private(Mat_SubSppt *submatj) 2414 { 2415 PetscErrorCode ierr; 2416 PetscInt i; 2417 2418 PetscFunctionBegin; 2419 if (!submatj->id) { /* delete data that are linked only to submats[id=0] */ 2420 ierr = PetscFree4(submatj->sbuf1,submatj->ptr,submatj->tmp,submatj->ctr);CHKERRQ(ierr); 2421 2422 for (i=0; i<submatj->nrqr; ++i) { 2423 ierr = PetscFree(submatj->sbuf2[i]);CHKERRQ(ierr); 2424 } 2425 ierr = PetscFree3(submatj->sbuf2,submatj->req_size,submatj->req_source1);CHKERRQ(ierr); 2426 2427 if (submatj->rbuf1) { 2428 ierr = PetscFree(submatj->rbuf1[0]);CHKERRQ(ierr); 2429 ierr = PetscFree(submatj->rbuf1);CHKERRQ(ierr); 2430 } 2431 2432 for (i=0; i<submatj->nrqs; ++i) { 2433 ierr = PetscFree(submatj->rbuf3[i]);CHKERRQ(ierr); 2434 } 2435 ierr = PetscFree3(submatj->req_source2,submatj->rbuf2,submatj->rbuf3);CHKERRQ(ierr); 2436 ierr = PetscFree(submatj->pa);CHKERRQ(ierr); 2437 } 2438 2439 #if defined(PETSC_USE_CTABLE) 2440 ierr = PetscTableDestroy((PetscTable*)&submatj->rmap);CHKERRQ(ierr); 2441 if (submatj->cmap_loc) {ierr = PetscFree(submatj->cmap_loc);CHKERRQ(ierr);} 2442 ierr = PetscFree(submatj->rmap_loc);CHKERRQ(ierr); 2443 #else 2444 ierr = PetscFree(submatj->rmap);CHKERRQ(ierr); 2445 #endif 2446 2447 if (!submatj->allcolumns) { 2448 #if defined(PETSC_USE_CTABLE) 2449 ierr = PetscTableDestroy((PetscTable*)&submatj->cmap);CHKERRQ(ierr); 2450 #else 2451 ierr = PetscFree(submatj->cmap);CHKERRQ(ierr); 2452 #endif 2453 } 2454 ierr = PetscFree(submatj->row2proc);CHKERRQ(ierr); 2455 2456 ierr = PetscFree(submatj);CHKERRQ(ierr); 2457 PetscFunctionReturn(0); 2458 } 2459 2460 PetscErrorCode MatDestroy_SeqAIJ_Submatrices(Mat C) 2461 { 2462 PetscErrorCode ierr; 2463 Mat_SeqAIJ *c = (Mat_SeqAIJ*)C->data; 2464 Mat_SubSppt *submatj = c->submatis1; 2465 2466 PetscFunctionBegin; 2467 ierr = submatj->destroy(C);CHKERRQ(ierr); 2468 ierr = MatDestroySubMatrices_Private(submatj);CHKERRQ(ierr); 2469 PetscFunctionReturn(0); 2470 } 2471 2472 PetscErrorCode MatCreateSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2473 { 2474 PetscErrorCode ierr; 2475 PetscInt i; 2476 2477 PetscFunctionBegin; 2478 if (scall == MAT_INITIAL_MATRIX) { 2479 ierr = PetscCalloc1(n+1,B);CHKERRQ(ierr); 2480 } 2481 2482 for (i=0; i<n; i++) { 2483 ierr = MatCreateSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2484 } 2485 PetscFunctionReturn(0); 2486 } 2487 2488 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2489 { 2490 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2491 PetscErrorCode ierr; 2492 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2493 const PetscInt *idx; 2494 PetscInt start,end,*ai,*aj; 2495 PetscBT table; 2496 2497 PetscFunctionBegin; 2498 m = A->rmap->n; 2499 ai = a->i; 2500 aj = a->j; 2501 2502 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2503 2504 ierr = PetscMalloc1(m+1,&nidx);CHKERRQ(ierr); 2505 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2506 2507 for (i=0; i<is_max; i++) { 2508 /* Initialize the two local arrays */ 2509 isz = 0; 2510 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2511 2512 /* Extract the indices, assume there can be duplicate entries */ 2513 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2514 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2515 2516 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2517 for (j=0; j<n; ++j) { 2518 if (!PetscBTLookupSet(table,idx[j])) nidx[isz++] = idx[j]; 2519 } 2520 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2521 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2522 2523 k = 0; 2524 for (j=0; j<ov; j++) { /* for each overlap */ 2525 n = isz; 2526 for (; k<n; k++) { /* do only those rows in nidx[k], which are not done yet */ 2527 row = nidx[k]; 2528 start = ai[row]; 2529 end = ai[row+1]; 2530 for (l = start; l<end; l++) { 2531 val = aj[l]; 2532 if (!PetscBTLookupSet(table,val)) nidx[isz++] = val; 2533 } 2534 } 2535 } 2536 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2537 } 2538 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2539 ierr = PetscFree(nidx);CHKERRQ(ierr); 2540 PetscFunctionReturn(0); 2541 } 2542 2543 /* -------------------------------------------------------------- */ 2544 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2545 { 2546 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2547 PetscErrorCode ierr; 2548 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2549 const PetscInt *row,*col; 2550 PetscInt *cnew,j,*lens; 2551 IS icolp,irowp; 2552 PetscInt *cwork = NULL; 2553 PetscScalar *vwork = NULL; 2554 2555 PetscFunctionBegin; 2556 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2557 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2558 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2559 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2560 2561 /* determine lengths of permuted rows */ 2562 ierr = PetscMalloc1(m+1,&lens);CHKERRQ(ierr); 2563 for (i=0; i<m; i++) lens[row[i]] = a->i[i+1] - a->i[i]; 2564 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 2565 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2566 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 2567 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2568 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2569 ierr = PetscFree(lens);CHKERRQ(ierr); 2570 2571 ierr = PetscMalloc1(n,&cnew);CHKERRQ(ierr); 2572 for (i=0; i<m; i++) { 2573 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2574 for (j=0; j<nz; j++) cnew[j] = col[cwork[j]]; 2575 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2576 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2577 } 2578 ierr = PetscFree(cnew);CHKERRQ(ierr); 2579 2580 (*B)->assembled = PETSC_FALSE; 2581 2582 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2583 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2584 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2585 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2586 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2587 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2588 PetscFunctionReturn(0); 2589 } 2590 2591 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2592 { 2593 PetscErrorCode ierr; 2594 2595 PetscFunctionBegin; 2596 /* If the two matrices have the same copy implementation, use fast copy. */ 2597 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2598 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2599 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2600 2601 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"); 2602 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2603 } else { 2604 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2605 } 2606 PetscFunctionReturn(0); 2607 } 2608 2609 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2610 { 2611 PetscErrorCode ierr; 2612 2613 PetscFunctionBegin; 2614 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2615 PetscFunctionReturn(0); 2616 } 2617 2618 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2619 { 2620 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2621 2622 PetscFunctionBegin; 2623 *array = a->a; 2624 PetscFunctionReturn(0); 2625 } 2626 2627 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2628 { 2629 PetscFunctionBegin; 2630 PetscFunctionReturn(0); 2631 } 2632 2633 /* 2634 Computes the number of nonzeros per row needed for preallocation when X and Y 2635 have different nonzero structure. 2636 */ 2637 PetscErrorCode MatAXPYGetPreallocation_SeqX_private(PetscInt m,const PetscInt *xi,const PetscInt *xj,const PetscInt *yi,const PetscInt *yj,PetscInt *nnz) 2638 { 2639 PetscInt i,j,k,nzx,nzy; 2640 2641 PetscFunctionBegin; 2642 /* Set the number of nonzeros in the new matrix */ 2643 for (i=0; i<m; i++) { 2644 const PetscInt *xjj = xj+xi[i],*yjj = yj+yi[i]; 2645 nzx = xi[i+1] - xi[i]; 2646 nzy = yi[i+1] - yi[i]; 2647 nnz[i] = 0; 2648 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2649 for (; k<nzy && yjj[k]<xjj[j]; k++) nnz[i]++; /* Catch up to X */ 2650 if (k<nzy && yjj[k]==xjj[j]) k++; /* Skip duplicate */ 2651 nnz[i]++; 2652 } 2653 for (; k<nzy; k++) nnz[i]++; 2654 } 2655 PetscFunctionReturn(0); 2656 } 2657 2658 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt *nnz) 2659 { 2660 PetscInt m = Y->rmap->N; 2661 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2662 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2663 PetscErrorCode ierr; 2664 2665 PetscFunctionBegin; 2666 /* Set the number of nonzeros in the new matrix */ 2667 ierr = MatAXPYGetPreallocation_SeqX_private(m,x->i,x->j,y->i,y->j,nnz);CHKERRQ(ierr); 2668 PetscFunctionReturn(0); 2669 } 2670 2671 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2672 { 2673 PetscErrorCode ierr; 2674 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data,*y = (Mat_SeqAIJ*)Y->data; 2675 PetscBLASInt one=1,bnz; 2676 2677 PetscFunctionBegin; 2678 ierr = PetscBLASIntCast(x->nz,&bnz);CHKERRQ(ierr); 2679 if (str == SAME_NONZERO_PATTERN) { 2680 PetscScalar alpha = a; 2681 PetscStackCallBLAS("BLASaxpy",BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one)); 2682 ierr = MatSeqAIJInvalidateDiagonal(Y);CHKERRQ(ierr); 2683 ierr = PetscObjectStateIncrease((PetscObject)Y);CHKERRQ(ierr); 2684 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2685 ierr = MatAXPY_Basic(Y,a,X,str);CHKERRQ(ierr); 2686 } else { 2687 Mat B; 2688 PetscInt *nnz; 2689 ierr = PetscMalloc1(Y->rmap->N,&nnz);CHKERRQ(ierr); 2690 ierr = MatCreate(PetscObjectComm((PetscObject)Y),&B);CHKERRQ(ierr); 2691 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2692 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2693 ierr = MatSetBlockSizesFromMats(B,Y,Y);CHKERRQ(ierr); 2694 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2695 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2696 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2697 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2698 ierr = MatHeaderReplace(Y,&B);CHKERRQ(ierr); 2699 ierr = PetscFree(nnz);CHKERRQ(ierr); 2700 } 2701 PetscFunctionReturn(0); 2702 } 2703 2704 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2705 { 2706 #if defined(PETSC_USE_COMPLEX) 2707 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 2708 PetscInt i,nz; 2709 PetscScalar *a; 2710 2711 PetscFunctionBegin; 2712 nz = aij->nz; 2713 a = aij->a; 2714 for (i=0; i<nz; i++) a[i] = PetscConj(a[i]); 2715 #else 2716 PetscFunctionBegin; 2717 #endif 2718 PetscFunctionReturn(0); 2719 } 2720 2721 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2722 { 2723 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2724 PetscErrorCode ierr; 2725 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2726 PetscReal atmp; 2727 PetscScalar *x; 2728 MatScalar *aa; 2729 2730 PetscFunctionBegin; 2731 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2732 aa = a->a; 2733 ai = a->i; 2734 aj = a->j; 2735 2736 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2737 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2738 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2739 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2740 for (i=0; i<m; i++) { 2741 ncols = ai[1] - ai[0]; ai++; 2742 x[i] = 0.0; 2743 for (j=0; j<ncols; j++) { 2744 atmp = PetscAbsScalar(*aa); 2745 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2746 aa++; aj++; 2747 } 2748 } 2749 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2750 PetscFunctionReturn(0); 2751 } 2752 2753 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2754 { 2755 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2756 PetscErrorCode ierr; 2757 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2758 PetscScalar *x; 2759 MatScalar *aa; 2760 2761 PetscFunctionBegin; 2762 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2763 aa = a->a; 2764 ai = a->i; 2765 aj = a->j; 2766 2767 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2768 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2769 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2770 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2771 for (i=0; i<m; i++) { 2772 ncols = ai[1] - ai[0]; ai++; 2773 if (ncols == A->cmap->n) { /* row is dense */ 2774 x[i] = *aa; if (idx) idx[i] = 0; 2775 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2776 x[i] = 0.0; 2777 if (idx) { 2778 idx[i] = 0; /* in case ncols is zero */ 2779 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2780 if (aj[j] > j) { 2781 idx[i] = j; 2782 break; 2783 } 2784 } 2785 } 2786 } 2787 for (j=0; j<ncols; j++) { 2788 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2789 aa++; aj++; 2790 } 2791 } 2792 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2793 PetscFunctionReturn(0); 2794 } 2795 2796 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2797 { 2798 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2799 PetscErrorCode ierr; 2800 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2801 PetscReal atmp; 2802 PetscScalar *x; 2803 MatScalar *aa; 2804 2805 PetscFunctionBegin; 2806 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2807 aa = a->a; 2808 ai = a->i; 2809 aj = a->j; 2810 2811 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2812 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2813 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2814 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); 2815 for (i=0; i<m; i++) { 2816 ncols = ai[1] - ai[0]; ai++; 2817 if (ncols) { 2818 /* Get first nonzero */ 2819 for (j = 0; j < ncols; j++) { 2820 atmp = PetscAbsScalar(aa[j]); 2821 if (atmp > 1.0e-12) { 2822 x[i] = atmp; 2823 if (idx) idx[i] = aj[j]; 2824 break; 2825 } 2826 } 2827 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2828 } else { 2829 x[i] = 0.0; if (idx) idx[i] = 0; 2830 } 2831 for (j = 0; j < ncols; j++) { 2832 atmp = PetscAbsScalar(*aa); 2833 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2834 aa++; aj++; 2835 } 2836 } 2837 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2838 PetscFunctionReturn(0); 2839 } 2840 2841 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2842 { 2843 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2844 PetscErrorCode ierr; 2845 PetscInt i,j,m = A->rmap->n,ncols,n; 2846 const PetscInt *ai,*aj; 2847 PetscScalar *x; 2848 const MatScalar *aa; 2849 2850 PetscFunctionBegin; 2851 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2852 aa = a->a; 2853 ai = a->i; 2854 aj = a->j; 2855 2856 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2857 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2858 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2859 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2860 for (i=0; i<m; i++) { 2861 ncols = ai[1] - ai[0]; ai++; 2862 if (ncols == A->cmap->n) { /* row is dense */ 2863 x[i] = *aa; if (idx) idx[i] = 0; 2864 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2865 x[i] = 0.0; 2866 if (idx) { /* find first implicit 0.0 in the row */ 2867 idx[i] = 0; /* in case ncols is zero */ 2868 for (j=0; j<ncols; j++) { 2869 if (aj[j] > j) { 2870 idx[i] = j; 2871 break; 2872 } 2873 } 2874 } 2875 } 2876 for (j=0; j<ncols; j++) { 2877 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2878 aa++; aj++; 2879 } 2880 } 2881 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2882 PetscFunctionReturn(0); 2883 } 2884 2885 #include <petscblaslapack.h> 2886 #include <petsc/private/kernels/blockinvert.h> 2887 2888 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2889 { 2890 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2891 PetscErrorCode ierr; 2892 PetscInt i,bs = PetscAbs(A->rmap->bs),mbs = A->rmap->n/bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2893 MatScalar *diag,work[25],*v_work; 2894 PetscReal shift = 0.0; 2895 PetscBool allowzeropivot,zeropivotdetected=PETSC_FALSE; 2896 2897 PetscFunctionBegin; 2898 allowzeropivot = PetscNot(A->erroriffailure); 2899 if (a->ibdiagvalid) { 2900 if (values) *values = a->ibdiag; 2901 PetscFunctionReturn(0); 2902 } 2903 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2904 if (!a->ibdiag) { 2905 ierr = PetscMalloc1(bs2*mbs,&a->ibdiag);CHKERRQ(ierr); 2906 ierr = PetscLogObjectMemory((PetscObject)A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2907 } 2908 diag = a->ibdiag; 2909 if (values) *values = a->ibdiag; 2910 /* factor and invert each block */ 2911 switch (bs) { 2912 case 1: 2913 for (i=0; i<mbs; i++) { 2914 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2915 if (PetscAbsScalar(diag[i] + shift) < PETSC_MACHINE_EPSILON) { 2916 if (allowzeropivot) { 2917 A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2918 A->factorerror_zeropivot_value = PetscAbsScalar(diag[i]); 2919 A->factorerror_zeropivot_row = i; 2920 ierr = PetscInfo3(A,"Zero pivot, row %D pivot %g tolerance %g\n",i,(double)PetscAbsScalar(diag[i]),(double)PETSC_MACHINE_EPSILON);CHKERRQ(ierr); 2921 } 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); 2922 } 2923 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2924 } 2925 break; 2926 case 2: 2927 for (i=0; i<mbs; i++) { 2928 ij[0] = 2*i; ij[1] = 2*i + 1; 2929 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2930 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2931 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2932 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2933 diag += 4; 2934 } 2935 break; 2936 case 3: 2937 for (i=0; i<mbs; i++) { 2938 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2939 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2940 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2941 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2942 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2943 diag += 9; 2944 } 2945 break; 2946 case 4: 2947 for (i=0; i<mbs; i++) { 2948 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2949 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2950 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2951 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2952 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2953 diag += 16; 2954 } 2955 break; 2956 case 5: 2957 for (i=0; i<mbs; i++) { 2958 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2959 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2960 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2961 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2962 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 2963 diag += 25; 2964 } 2965 break; 2966 case 6: 2967 for (i=0; i<mbs; i++) { 2968 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; 2969 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2970 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2971 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2972 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 2973 diag += 36; 2974 } 2975 break; 2976 case 7: 2977 for (i=0; i<mbs; i++) { 2978 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; 2979 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 2980 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2981 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2982 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 2983 diag += 49; 2984 } 2985 break; 2986 default: 2987 ierr = PetscMalloc3(bs,&v_work,bs,&v_pivots,bs,&IJ);CHKERRQ(ierr); 2988 for (i=0; i<mbs; i++) { 2989 for (j=0; j<bs; j++) { 2990 IJ[j] = bs*i + j; 2991 } 2992 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 2993 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work,allowzeropivot,&zeropivotdetected);CHKERRQ(ierr); 2994 if (zeropivotdetected) A->factorerrortype = MAT_FACTOR_NUMERIC_ZEROPIVOT; 2995 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 2996 diag += bs2; 2997 } 2998 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 2999 } 3000 a->ibdiagvalid = PETSC_TRUE; 3001 PetscFunctionReturn(0); 3002 } 3003 3004 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3005 { 3006 PetscErrorCode ierr; 3007 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3008 PetscScalar a; 3009 PetscInt m,n,i,j,col; 3010 3011 PetscFunctionBegin; 3012 if (!x->assembled) { 3013 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3014 for (i=0; i<m; i++) { 3015 for (j=0; j<aij->imax[i]; j++) { 3016 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3017 col = (PetscInt)(n*PetscRealPart(a)); 3018 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3019 } 3020 } 3021 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3022 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3023 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3024 PetscFunctionReturn(0); 3025 } 3026 3027 PetscErrorCode MatShift_SeqAIJ(Mat Y,PetscScalar a) 3028 { 3029 PetscErrorCode ierr; 3030 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)Y->data; 3031 3032 PetscFunctionBegin; 3033 if (!Y->preallocated || !aij->nz) { 3034 ierr = MatSeqAIJSetPreallocation(Y,1,NULL);CHKERRQ(ierr); 3035 } 3036 ierr = MatShift_Basic(Y,a);CHKERRQ(ierr); 3037 PetscFunctionReturn(0); 3038 } 3039 3040 /* -------------------------------------------------------------------*/ 3041 static struct _MatOps MatOps_Values = { MatSetValues_SeqAIJ, 3042 MatGetRow_SeqAIJ, 3043 MatRestoreRow_SeqAIJ, 3044 MatMult_SeqAIJ, 3045 /* 4*/ MatMultAdd_SeqAIJ, 3046 MatMultTranspose_SeqAIJ, 3047 MatMultTransposeAdd_SeqAIJ, 3048 0, 3049 0, 3050 0, 3051 /* 10*/ 0, 3052 MatLUFactor_SeqAIJ, 3053 0, 3054 MatSOR_SeqAIJ, 3055 MatTranspose_SeqAIJ, 3056 /*1 5*/ MatGetInfo_SeqAIJ, 3057 MatEqual_SeqAIJ, 3058 MatGetDiagonal_SeqAIJ, 3059 MatDiagonalScale_SeqAIJ, 3060 MatNorm_SeqAIJ, 3061 /* 20*/ 0, 3062 MatAssemblyEnd_SeqAIJ, 3063 MatSetOption_SeqAIJ, 3064 MatZeroEntries_SeqAIJ, 3065 /* 24*/ MatZeroRows_SeqAIJ, 3066 0, 3067 0, 3068 0, 3069 0, 3070 /* 29*/ MatSetUp_SeqAIJ, 3071 0, 3072 0, 3073 0, 3074 0, 3075 /* 34*/ MatDuplicate_SeqAIJ, 3076 0, 3077 0, 3078 MatILUFactor_SeqAIJ, 3079 0, 3080 /* 39*/ MatAXPY_SeqAIJ, 3081 MatCreateSubMatrices_SeqAIJ, 3082 MatIncreaseOverlap_SeqAIJ, 3083 MatGetValues_SeqAIJ, 3084 MatCopy_SeqAIJ, 3085 /* 44*/ MatGetRowMax_SeqAIJ, 3086 MatScale_SeqAIJ, 3087 MatShift_SeqAIJ, 3088 MatDiagonalSet_SeqAIJ, 3089 MatZeroRowsColumns_SeqAIJ, 3090 /* 49*/ MatSetRandom_SeqAIJ, 3091 MatGetRowIJ_SeqAIJ, 3092 MatRestoreRowIJ_SeqAIJ, 3093 MatGetColumnIJ_SeqAIJ, 3094 MatRestoreColumnIJ_SeqAIJ, 3095 /* 54*/ MatFDColoringCreate_SeqXAIJ, 3096 0, 3097 0, 3098 MatPermute_SeqAIJ, 3099 0, 3100 /* 59*/ 0, 3101 MatDestroy_SeqAIJ, 3102 MatView_SeqAIJ, 3103 0, 3104 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3105 /* 64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3106 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3107 0, 3108 0, 3109 0, 3110 /* 69*/ MatGetRowMaxAbs_SeqAIJ, 3111 MatGetRowMinAbs_SeqAIJ, 3112 0, 3113 0, 3114 0, 3115 /* 74*/ 0, 3116 MatFDColoringApply_AIJ, 3117 0, 3118 0, 3119 0, 3120 /* 79*/ MatFindZeroDiagonals_SeqAIJ, 3121 0, 3122 0, 3123 0, 3124 MatLoad_SeqAIJ, 3125 /* 84*/ MatIsSymmetric_SeqAIJ, 3126 MatIsHermitian_SeqAIJ, 3127 0, 3128 0, 3129 0, 3130 /* 89*/ MatMatMult_SeqAIJ_SeqAIJ, 3131 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3132 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3133 MatPtAP_SeqAIJ_SeqAIJ, 3134 MatPtAPSymbolic_SeqAIJ_SeqAIJ_DenseAxpy, 3135 /* 94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3136 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3137 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3138 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3139 0, 3140 /* 99*/ 0, 3141 0, 3142 0, 3143 MatConjugate_SeqAIJ, 3144 0, 3145 /*104*/ MatSetValuesRow_SeqAIJ, 3146 MatRealPart_SeqAIJ, 3147 MatImaginaryPart_SeqAIJ, 3148 0, 3149 0, 3150 /*109*/ MatMatSolve_SeqAIJ, 3151 0, 3152 MatGetRowMin_SeqAIJ, 3153 0, 3154 MatMissingDiagonal_SeqAIJ, 3155 /*114*/ 0, 3156 0, 3157 0, 3158 0, 3159 0, 3160 /*119*/ 0, 3161 0, 3162 0, 3163 0, 3164 MatGetMultiProcBlock_SeqAIJ, 3165 /*124*/ MatFindNonzeroRows_SeqAIJ, 3166 MatGetColumnNorms_SeqAIJ, 3167 MatInvertBlockDiagonal_SeqAIJ, 3168 0, 3169 0, 3170 /*129*/ 0, 3171 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3172 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3173 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3174 MatTransposeColoringCreate_SeqAIJ, 3175 /*134*/ MatTransColoringApplySpToDen_SeqAIJ, 3176 MatTransColoringApplyDenToSp_SeqAIJ, 3177 MatRARt_SeqAIJ_SeqAIJ, 3178 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3179 MatRARtNumeric_SeqAIJ_SeqAIJ, 3180 /*139*/0, 3181 0, 3182 0, 3183 MatFDColoringSetUp_SeqXAIJ, 3184 MatFindOffBlockDiagonalEntries_SeqAIJ, 3185 /*144*/MatCreateMPIMatConcatenateSeqMat_SeqAIJ 3186 }; 3187 3188 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3189 { 3190 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3191 PetscInt i,nz,n; 3192 3193 PetscFunctionBegin; 3194 nz = aij->maxnz; 3195 n = mat->rmap->n; 3196 for (i=0; i<nz; i++) { 3197 aij->j[i] = indices[i]; 3198 } 3199 aij->nz = nz; 3200 for (i=0; i<n; i++) { 3201 aij->ilen[i] = aij->imax[i]; 3202 } 3203 PetscFunctionReturn(0); 3204 } 3205 3206 /*@ 3207 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3208 in the matrix. 3209 3210 Input Parameters: 3211 + mat - the SeqAIJ matrix 3212 - indices - the column indices 3213 3214 Level: advanced 3215 3216 Notes: 3217 This can be called if you have precomputed the nonzero structure of the 3218 matrix and want to provide it to the matrix object to improve the performance 3219 of the MatSetValues() operation. 3220 3221 You MUST have set the correct numbers of nonzeros per row in the call to 3222 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3223 3224 MUST be called before any calls to MatSetValues(); 3225 3226 The indices should start with zero, not one. 3227 3228 @*/ 3229 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3230 { 3231 PetscErrorCode ierr; 3232 3233 PetscFunctionBegin; 3234 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3235 PetscValidPointer(indices,2); 3236 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt*),(mat,indices));CHKERRQ(ierr); 3237 PetscFunctionReturn(0); 3238 } 3239 3240 /* ----------------------------------------------------------------------------------------*/ 3241 3242 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3243 { 3244 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3245 PetscErrorCode ierr; 3246 size_t nz = aij->i[mat->rmap->n]; 3247 3248 PetscFunctionBegin; 3249 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3250 3251 /* allocate space for values if not already there */ 3252 if (!aij->saved_values) { 3253 ierr = PetscMalloc1(nz+1,&aij->saved_values);CHKERRQ(ierr); 3254 ierr = PetscLogObjectMemory((PetscObject)mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3255 } 3256 3257 /* copy values over */ 3258 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3259 PetscFunctionReturn(0); 3260 } 3261 3262 /*@ 3263 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3264 example, reuse of the linear part of a Jacobian, while recomputing the 3265 nonlinear portion. 3266 3267 Collect on Mat 3268 3269 Input Parameters: 3270 . mat - the matrix (currently only AIJ matrices support this option) 3271 3272 Level: advanced 3273 3274 Common Usage, with SNESSolve(): 3275 $ Create Jacobian matrix 3276 $ Set linear terms into matrix 3277 $ Apply boundary conditions to matrix, at this time matrix must have 3278 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3279 $ boundary conditions again will not change the nonzero structure 3280 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3281 $ ierr = MatStoreValues(mat); 3282 $ Call SNESSetJacobian() with matrix 3283 $ In your Jacobian routine 3284 $ ierr = MatRetrieveValues(mat); 3285 $ Set nonlinear terms in matrix 3286 3287 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3288 $ // build linear portion of Jacobian 3289 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3290 $ ierr = MatStoreValues(mat); 3291 $ loop over nonlinear iterations 3292 $ ierr = MatRetrieveValues(mat); 3293 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3294 $ // call MatAssemblyBegin/End() on matrix 3295 $ Solve linear system with Jacobian 3296 $ endloop 3297 3298 Notes: 3299 Matrix must already be assemblied before calling this routine 3300 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3301 calling this routine. 3302 3303 When this is called multiple times it overwrites the previous set of stored values 3304 and does not allocated additional space. 3305 3306 .seealso: MatRetrieveValues() 3307 3308 @*/ 3309 PetscErrorCode MatStoreValues(Mat mat) 3310 { 3311 PetscErrorCode ierr; 3312 3313 PetscFunctionBegin; 3314 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3315 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3316 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3317 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3318 PetscFunctionReturn(0); 3319 } 3320 3321 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3322 { 3323 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)mat->data; 3324 PetscErrorCode ierr; 3325 PetscInt nz = aij->i[mat->rmap->n]; 3326 3327 PetscFunctionBegin; 3328 if (!aij->nonew) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3329 if (!aij->saved_values) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3330 /* copy values over */ 3331 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3332 PetscFunctionReturn(0); 3333 } 3334 3335 /*@ 3336 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3337 example, reuse of the linear part of a Jacobian, while recomputing the 3338 nonlinear portion. 3339 3340 Collect on Mat 3341 3342 Input Parameters: 3343 . mat - the matrix (currently only AIJ matrices support this option) 3344 3345 Level: advanced 3346 3347 .seealso: MatStoreValues() 3348 3349 @*/ 3350 PetscErrorCode MatRetrieveValues(Mat mat) 3351 { 3352 PetscErrorCode ierr; 3353 3354 PetscFunctionBegin; 3355 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3356 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3357 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3358 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3359 PetscFunctionReturn(0); 3360 } 3361 3362 3363 /* --------------------------------------------------------------------------------*/ 3364 /*@C 3365 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3366 (the default parallel PETSc format). For good matrix assembly performance 3367 the user should preallocate the matrix storage by setting the parameter nz 3368 (or the array nnz). By setting these parameters accurately, performance 3369 during matrix assembly can be increased by more than a factor of 50. 3370 3371 Collective on MPI_Comm 3372 3373 Input Parameters: 3374 + comm - MPI communicator, set to PETSC_COMM_SELF 3375 . m - number of rows 3376 . n - number of columns 3377 . nz - number of nonzeros per row (same for all rows) 3378 - nnz - array containing the number of nonzeros in the various rows 3379 (possibly different for each row) or NULL 3380 3381 Output Parameter: 3382 . A - the matrix 3383 3384 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3385 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3386 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3387 3388 Notes: 3389 If nnz is given then nz is ignored 3390 3391 The AIJ format (also called the Yale sparse matrix format or 3392 compressed row storage), is fully compatible with standard Fortran 77 3393 storage. That is, the stored row and column indices can begin at 3394 either one (as in Fortran) or zero. See the users' manual for details. 3395 3396 Specify the preallocated storage with either nz or nnz (not both). 3397 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3398 allocation. For large problems you MUST preallocate memory or you 3399 will get TERRIBLE performance, see the users' manual chapter on matrices. 3400 3401 By default, this format uses inodes (identical nodes) when possible, to 3402 improve numerical efficiency of matrix-vector products and solves. We 3403 search for consecutive rows with the same nonzero structure, thereby 3404 reusing matrix information to achieve increased efficiency. 3405 3406 Options Database Keys: 3407 + -mat_no_inode - Do not use inodes 3408 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3409 3410 Level: intermediate 3411 3412 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3413 3414 @*/ 3415 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3416 { 3417 PetscErrorCode ierr; 3418 3419 PetscFunctionBegin; 3420 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3421 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3422 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3423 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3424 PetscFunctionReturn(0); 3425 } 3426 3427 /*@C 3428 MatSeqAIJSetPreallocation - For good matrix assembly performance 3429 the user should preallocate the matrix storage by setting the parameter nz 3430 (or the array nnz). By setting these parameters accurately, performance 3431 during matrix assembly can be increased by more than a factor of 50. 3432 3433 Collective on MPI_Comm 3434 3435 Input Parameters: 3436 + B - The matrix 3437 . nz - number of nonzeros per row (same for all rows) 3438 - nnz - array containing the number of nonzeros in the various rows 3439 (possibly different for each row) or NULL 3440 3441 Notes: 3442 If nnz is given then nz is ignored 3443 3444 The AIJ format (also called the Yale sparse matrix format or 3445 compressed row storage), is fully compatible with standard Fortran 77 3446 storage. That is, the stored row and column indices can begin at 3447 either one (as in Fortran) or zero. See the users' manual for details. 3448 3449 Specify the preallocated storage with either nz or nnz (not both). 3450 Set nz=PETSC_DEFAULT and nnz=NULL for PETSc to control dynamic memory 3451 allocation. For large problems you MUST preallocate memory or you 3452 will get TERRIBLE performance, see the users' manual chapter on matrices. 3453 3454 You can call MatGetInfo() to get information on how effective the preallocation was; 3455 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3456 You can also run with the option -info and look for messages with the string 3457 malloc in them to see if additional memory allocation was needed. 3458 3459 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3460 entries or columns indices 3461 3462 By default, this format uses inodes (identical nodes) when possible, to 3463 improve numerical efficiency of matrix-vector products and solves. We 3464 search for consecutive rows with the same nonzero structure, thereby 3465 reusing matrix information to achieve increased efficiency. 3466 3467 Options Database Keys: 3468 + -mat_no_inode - Do not use inodes 3469 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3470 - -mat_aij_oneindex - Internally use indexing starting at 1 3471 rather than 0. Note that when calling MatSetValues(), 3472 the user still MUST index entries starting at 0! 3473 3474 Level: intermediate 3475 3476 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3477 3478 @*/ 3479 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3480 { 3481 PetscErrorCode ierr; 3482 3483 PetscFunctionBegin; 3484 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3485 PetscValidType(B,1); 3486 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3487 PetscFunctionReturn(0); 3488 } 3489 3490 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3491 { 3492 Mat_SeqAIJ *b; 3493 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3494 PetscErrorCode ierr; 3495 PetscInt i; 3496 3497 PetscFunctionBegin; 3498 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3499 if (nz == MAT_SKIP_ALLOCATION) { 3500 skipallocation = PETSC_TRUE; 3501 nz = 0; 3502 } 3503 3504 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3505 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3506 3507 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3508 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %D",nz); 3509 if (nnz) { 3510 for (i=0; i<B->rmap->n; i++) { 3511 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]); 3512 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); 3513 } 3514 } 3515 3516 B->preallocated = PETSC_TRUE; 3517 3518 b = (Mat_SeqAIJ*)B->data; 3519 3520 if (!skipallocation) { 3521 if (!b->imax) { 3522 ierr = PetscMalloc2(B->rmap->n,&b->imax,B->rmap->n,&b->ilen);CHKERRQ(ierr); 3523 ierr = PetscLogObjectMemory((PetscObject)B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3524 } 3525 if (!nnz) { 3526 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3527 else if (nz < 0) nz = 1; 3528 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3529 nz = nz*B->rmap->n; 3530 } else { 3531 nz = 0; 3532 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3533 } 3534 /* b->ilen will count nonzeros in each row so far. */ 3535 for (i=0; i<B->rmap->n; i++) b->ilen[i] = 0; 3536 3537 /* allocate the matrix space */ 3538 /* FIXME: should B's old memory be unlogged? */ 3539 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3540 ierr = PetscMalloc3(nz,&b->a,nz,&b->j,B->rmap->n+1,&b->i);CHKERRQ(ierr); 3541 ierr = PetscLogObjectMemory((PetscObject)B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3542 b->i[0] = 0; 3543 for (i=1; i<B->rmap->n+1; i++) { 3544 b->i[i] = b->i[i-1] + b->imax[i-1]; 3545 } 3546 b->singlemalloc = PETSC_TRUE; 3547 b->free_a = PETSC_TRUE; 3548 b->free_ij = PETSC_TRUE; 3549 } else { 3550 b->free_a = PETSC_FALSE; 3551 b->free_ij = PETSC_FALSE; 3552 } 3553 3554 b->nz = 0; 3555 b->maxnz = nz; 3556 B->info.nz_unneeded = (double)b->maxnz; 3557 if (realalloc) { 3558 ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3559 } 3560 B->was_assembled = PETSC_FALSE; 3561 B->assembled = PETSC_FALSE; 3562 PetscFunctionReturn(0); 3563 } 3564 3565 /*@ 3566 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3567 3568 Input Parameters: 3569 + B - the matrix 3570 . i - the indices into j for the start of each row (starts with zero) 3571 . j - the column indices for each row (starts with zero) these must be sorted for each row 3572 - v - optional values in the matrix 3573 3574 Level: developer 3575 3576 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3577 3578 .keywords: matrix, aij, compressed row, sparse, sequential 3579 3580 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3581 @*/ 3582 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3583 { 3584 PetscErrorCode ierr; 3585 3586 PetscFunctionBegin; 3587 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3588 PetscValidType(B,1); 3589 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3590 PetscFunctionReturn(0); 3591 } 3592 3593 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3594 { 3595 PetscInt i; 3596 PetscInt m,n; 3597 PetscInt nz; 3598 PetscInt *nnz, nz_max = 0; 3599 PetscScalar *values; 3600 PetscErrorCode ierr; 3601 3602 PetscFunctionBegin; 3603 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3604 3605 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3606 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3607 3608 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3609 ierr = PetscMalloc1(m+1, &nnz);CHKERRQ(ierr); 3610 for (i = 0; i < m; i++) { 3611 nz = Ii[i+1]- Ii[i]; 3612 nz_max = PetscMax(nz_max, nz); 3613 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3614 nnz[i] = nz; 3615 } 3616 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3617 ierr = PetscFree(nnz);CHKERRQ(ierr); 3618 3619 if (v) { 3620 values = (PetscScalar*) v; 3621 } else { 3622 ierr = PetscCalloc1(nz_max, &values);CHKERRQ(ierr); 3623 } 3624 3625 for (i = 0; i < m; i++) { 3626 nz = Ii[i+1] - Ii[i]; 3627 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3628 } 3629 3630 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3631 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3632 3633 if (!v) { 3634 ierr = PetscFree(values);CHKERRQ(ierr); 3635 } 3636 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3637 PetscFunctionReturn(0); 3638 } 3639 3640 #include <../src/mat/impls/dense/seq/dense.h> 3641 #include <petsc/private/kernels/petscaxpy.h> 3642 3643 /* 3644 Computes (B'*A')' since computing B*A directly is untenable 3645 3646 n p p 3647 ( ) ( ) ( ) 3648 m ( A ) * n ( B ) = m ( C ) 3649 ( ) ( ) ( ) 3650 3651 */ 3652 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3653 { 3654 PetscErrorCode ierr; 3655 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3656 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3657 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3658 PetscInt i,n,m,q,p; 3659 const PetscInt *ii,*idx; 3660 const PetscScalar *b,*a,*a_q; 3661 PetscScalar *c,*c_q; 3662 3663 PetscFunctionBegin; 3664 m = A->rmap->n; 3665 n = A->cmap->n; 3666 p = B->cmap->n; 3667 a = sub_a->v; 3668 b = sub_b->a; 3669 c = sub_c->v; 3670 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3671 3672 ii = sub_b->i; 3673 idx = sub_b->j; 3674 for (i=0; i<n; i++) { 3675 q = ii[i+1] - ii[i]; 3676 while (q-->0) { 3677 c_q = c + m*(*idx); 3678 a_q = a + m*i; 3679 PetscKernelAXPY(c_q,*b,a_q,m); 3680 idx++; 3681 b++; 3682 } 3683 } 3684 PetscFunctionReturn(0); 3685 } 3686 3687 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3688 { 3689 PetscErrorCode ierr; 3690 PetscInt m=A->rmap->n,n=B->cmap->n; 3691 Mat Cmat; 3692 3693 PetscFunctionBegin; 3694 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); 3695 ierr = MatCreate(PetscObjectComm((PetscObject)A),&Cmat);CHKERRQ(ierr); 3696 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3697 ierr = MatSetBlockSizesFromMats(Cmat,A,B);CHKERRQ(ierr); 3698 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3699 ierr = MatSeqDenseSetPreallocation(Cmat,NULL);CHKERRQ(ierr); 3700 3701 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3702 3703 *C = Cmat; 3704 PetscFunctionReturn(0); 3705 } 3706 3707 /* ----------------------------------------------------------------*/ 3708 PETSC_INTERN PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3709 { 3710 PetscErrorCode ierr; 3711 3712 PetscFunctionBegin; 3713 if (scall == MAT_INITIAL_MATRIX) { 3714 ierr = PetscLogEventBegin(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3715 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3716 ierr = PetscLogEventEnd(MAT_MatMultSymbolic,A,B,0,0);CHKERRQ(ierr); 3717 } 3718 ierr = PetscLogEventBegin(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3719 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3720 ierr = PetscLogEventEnd(MAT_MatMultNumeric,A,B,0,0);CHKERRQ(ierr); 3721 PetscFunctionReturn(0); 3722 } 3723 3724 3725 /*MC 3726 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3727 based on compressed sparse row format. 3728 3729 Options Database Keys: 3730 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3731 3732 Level: beginner 3733 3734 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3735 M*/ 3736 3737 /*MC 3738 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3739 3740 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3741 and MATMPIAIJ otherwise. As a result, for single process communicators, 3742 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3743 for communicators controlling multiple processes. It is recommended that you call both of 3744 the above preallocation routines for simplicity. 3745 3746 Options Database Keys: 3747 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3748 3749 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3750 enough exist. 3751 3752 Level: beginner 3753 3754 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3755 M*/ 3756 3757 /*MC 3758 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3759 3760 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3761 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3762 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3763 for communicators controlling multiple processes. It is recommended that you call both of 3764 the above preallocation routines for simplicity. 3765 3766 Options Database Keys: 3767 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3768 3769 Level: beginner 3770 3771 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3772 M*/ 3773 3774 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3775 #if defined(PETSC_HAVE_ELEMENTAL) 3776 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_Elemental(Mat,MatType,MatReuse,Mat*); 3777 #endif 3778 #if defined(PETSC_HAVE_HYPRE) 3779 PETSC_INTERN PetscErrorCode MatConvert_AIJ_HYPRE(Mat A,MatType,MatReuse,Mat*); 3780 PETSC_INTERN PetscErrorCode MatMatMatMult_Transpose_AIJ_AIJ(Mat,Mat,Mat,MatReuse,PetscReal,Mat*); 3781 #endif 3782 PETSC_INTERN PetscErrorCode MatConvert_SeqAIJ_SeqDense(Mat,MatType,MatReuse,Mat*); 3783 3784 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3785 PETSC_EXTERN PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3786 PETSC_EXTERN PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3787 #endif 3788 3789 3790 /*@C 3791 MatSeqAIJGetArray - gives access to the array where the data for a MATSEQAIJ matrix is stored 3792 3793 Not Collective 3794 3795 Input Parameter: 3796 . mat - a MATSEQAIJ matrix 3797 3798 Output Parameter: 3799 . array - pointer to the data 3800 3801 Level: intermediate 3802 3803 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3804 @*/ 3805 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3806 { 3807 PetscErrorCode ierr; 3808 3809 PetscFunctionBegin; 3810 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3811 PetscFunctionReturn(0); 3812 } 3813 3814 /*@C 3815 MatSeqAIJGetMaxRowNonzeros - returns the maximum number of nonzeros in any row 3816 3817 Not Collective 3818 3819 Input Parameter: 3820 . mat - a MATSEQAIJ matrix 3821 3822 Output Parameter: 3823 . nz - the maximum number of nonzeros in any row 3824 3825 Level: intermediate 3826 3827 .seealso: MatSeqAIJRestoreArray(), MatSeqAIJGetArrayF90() 3828 @*/ 3829 PetscErrorCode MatSeqAIJGetMaxRowNonzeros(Mat A,PetscInt *nz) 3830 { 3831 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)A->data; 3832 3833 PetscFunctionBegin; 3834 *nz = aij->rmax; 3835 PetscFunctionReturn(0); 3836 } 3837 3838 /*@C 3839 MatSeqAIJRestoreArray - returns access to the array where the data for a MATSEQAIJ matrix is stored obtained by MatSeqAIJGetArray() 3840 3841 Not Collective 3842 3843 Input Parameters: 3844 . mat - a MATSEQAIJ matrix 3845 . array - pointer to the data 3846 3847 Level: intermediate 3848 3849 .seealso: MatSeqAIJGetArray(), MatSeqAIJRestoreArrayF90() 3850 @*/ 3851 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3852 { 3853 PetscErrorCode ierr; 3854 3855 PetscFunctionBegin; 3856 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3857 PetscFunctionReturn(0); 3858 } 3859 3860 PETSC_EXTERN PetscErrorCode MatCreate_SeqAIJ(Mat B) 3861 { 3862 Mat_SeqAIJ *b; 3863 PetscErrorCode ierr; 3864 PetscMPIInt size; 3865 3866 PetscFunctionBegin; 3867 ierr = MPI_Comm_size(PetscObjectComm((PetscObject)B),&size);CHKERRQ(ierr); 3868 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3869 3870 ierr = PetscNewLog(B,&b);CHKERRQ(ierr); 3871 3872 B->data = (void*)b; 3873 3874 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3875 3876 b->row = 0; 3877 b->col = 0; 3878 b->icol = 0; 3879 b->reallocs = 0; 3880 b->ignorezeroentries = PETSC_FALSE; 3881 b->roworiented = PETSC_TRUE; 3882 b->nonew = 0; 3883 b->diag = 0; 3884 b->solve_work = 0; 3885 B->spptr = 0; 3886 b->saved_values = 0; 3887 b->idiag = 0; 3888 b->mdiag = 0; 3889 b->ssor_work = 0; 3890 b->omega = 1.0; 3891 b->fshift = 0.0; 3892 b->idiagvalid = PETSC_FALSE; 3893 b->ibdiagvalid = PETSC_FALSE; 3894 b->keepnonzeropattern = PETSC_FALSE; 3895 3896 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3897 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJGetArray_C",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3898 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJRestoreArray_C",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3899 3900 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3901 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEnginePut_C",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3902 ierr = PetscObjectComposeFunction((PetscObject)B,"PetscMatlabEngineGet_C",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3903 #endif 3904 3905 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetColumnIndices_C",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3906 ierr = PetscObjectComposeFunction((PetscObject)B,"MatStoreValues_C",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3907 ierr = PetscObjectComposeFunction((PetscObject)B,"MatRetrieveValues_C",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3908 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqsbaij_C",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3909 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqbaij_C",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3910 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijperm_C",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3911 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3912 #if defined(PETSC_HAVE_ELEMENTAL) 3913 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_elemental_C",MatConvert_SeqAIJ_Elemental);CHKERRQ(ierr); 3914 #endif 3915 #if defined(PETSC_HAVE_HYPRE) 3916 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_hypre_C",MatConvert_AIJ_HYPRE);CHKERRQ(ierr); 3917 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMatMult_transpose_seqaij_seqaij_C",MatMatMatMult_Transpose_AIJ_AIJ);CHKERRQ(ierr); 3918 #endif 3919 ierr = PetscObjectComposeFunction((PetscObject)B,"MatConvert_seqaij_seqdense_C",MatConvert_SeqAIJ_SeqDense);CHKERRQ(ierr); 3920 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3921 ierr = PetscObjectComposeFunction((PetscObject)B,"MatIsHermitianTranspose_C",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3922 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocation_C",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3923 ierr = PetscObjectComposeFunction((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3924 ierr = PetscObjectComposeFunction((PetscObject)B,"MatReorderForNonzeroDiagonal_C",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3925 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMult_seqdense_seqaij_C",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3926 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3927 ierr = PetscObjectComposeFunction((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3928 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3929 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3930 PetscFunctionReturn(0); 3931 } 3932 3933 /* 3934 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 3935 */ 3936 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 3937 { 3938 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 3939 PetscErrorCode ierr; 3940 PetscInt i,m = A->rmap->n; 3941 3942 PetscFunctionBegin; 3943 c = (Mat_SeqAIJ*)C->data; 3944 3945 C->factortype = A->factortype; 3946 c->row = 0; 3947 c->col = 0; 3948 c->icol = 0; 3949 c->reallocs = 0; 3950 3951 C->assembled = PETSC_TRUE; 3952 3953 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 3954 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 3955 3956 ierr = PetscMalloc2(m,&c->imax,m,&c->ilen);CHKERRQ(ierr); 3957 ierr = PetscLogObjectMemory((PetscObject)C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 3958 for (i=0; i<m; i++) { 3959 c->imax[i] = a->imax[i]; 3960 c->ilen[i] = a->ilen[i]; 3961 } 3962 3963 /* allocate the matrix space */ 3964 if (mallocmatspace) { 3965 ierr = PetscMalloc3(a->i[m],&c->a,a->i[m],&c->j,m+1,&c->i);CHKERRQ(ierr); 3966 ierr = PetscLogObjectMemory((PetscObject)C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3967 3968 c->singlemalloc = PETSC_TRUE; 3969 3970 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3971 if (m > 0) { 3972 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 3973 if (cpvalues == MAT_COPY_VALUES) { 3974 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3975 } else { 3976 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 3977 } 3978 } 3979 } 3980 3981 c->ignorezeroentries = a->ignorezeroentries; 3982 c->roworiented = a->roworiented; 3983 c->nonew = a->nonew; 3984 if (a->diag) { 3985 ierr = PetscMalloc1(m+1,&c->diag);CHKERRQ(ierr); 3986 ierr = PetscLogObjectMemory((PetscObject)C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 3987 for (i=0; i<m; i++) { 3988 c->diag[i] = a->diag[i]; 3989 } 3990 } else c->diag = 0; 3991 3992 c->solve_work = 0; 3993 c->saved_values = 0; 3994 c->idiag = 0; 3995 c->ssor_work = 0; 3996 c->keepnonzeropattern = a->keepnonzeropattern; 3997 c->free_a = PETSC_TRUE; 3998 c->free_ij = PETSC_TRUE; 3999 4000 c->rmax = a->rmax; 4001 c->nz = a->nz; 4002 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4003 C->preallocated = PETSC_TRUE; 4004 4005 c->compressedrow.use = a->compressedrow.use; 4006 c->compressedrow.nrows = a->compressedrow.nrows; 4007 if (a->compressedrow.use) { 4008 i = a->compressedrow.nrows; 4009 ierr = PetscMalloc2(i+1,&c->compressedrow.i,i,&c->compressedrow.rindex);CHKERRQ(ierr); 4010 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4011 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4012 } else { 4013 c->compressedrow.use = PETSC_FALSE; 4014 c->compressedrow.i = NULL; 4015 c->compressedrow.rindex = NULL; 4016 } 4017 c->nonzerorowcnt = a->nonzerorowcnt; 4018 C->nonzerostate = A->nonzerostate; 4019 4020 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4021 ierr = PetscFunctionListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4022 PetscFunctionReturn(0); 4023 } 4024 4025 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4026 { 4027 PetscErrorCode ierr; 4028 4029 PetscFunctionBegin; 4030 ierr = MatCreate(PetscObjectComm((PetscObject)A),B);CHKERRQ(ierr); 4031 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4032 if (!(A->rmap->n % A->rmap->bs) && !(A->cmap->n % A->cmap->bs)) { 4033 ierr = MatSetBlockSizesFromMats(*B,A,A);CHKERRQ(ierr); 4034 } 4035 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4036 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4037 PetscFunctionReturn(0); 4038 } 4039 4040 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4041 { 4042 Mat_SeqAIJ *a; 4043 PetscErrorCode ierr; 4044 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4045 int fd; 4046 PetscMPIInt size; 4047 MPI_Comm comm; 4048 PetscInt bs = newMat->rmap->bs; 4049 4050 PetscFunctionBegin; 4051 /* force binary viewer to load .info file if it has not yet done so */ 4052 ierr = PetscViewerSetUp(viewer);CHKERRQ(ierr); 4053 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4054 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4055 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4056 4057 ierr = PetscOptionsBegin(comm,NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4058 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,NULL);CHKERRQ(ierr); 4059 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4060 if (bs < 0) bs = 1; 4061 ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr); 4062 4063 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4064 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4065 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4066 M = header[1]; N = header[2]; nz = header[3]; 4067 4068 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4069 4070 /* read in row lengths */ 4071 ierr = PetscMalloc1(M,&rowlengths);CHKERRQ(ierr); 4072 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4073 4074 /* check if sum of rowlengths is same as nz */ 4075 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4076 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); 4077 4078 /* set global size if not set already*/ 4079 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4080 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4081 } else { 4082 /* if sizes and type are already set, check if the matrix global sizes are correct */ 4083 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4084 if (rows < 0 && cols < 0) { /* user might provide local size instead of global size */ 4085 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4086 } 4087 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); 4088 } 4089 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4090 a = (Mat_SeqAIJ*)newMat->data; 4091 4092 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4093 4094 /* read in nonzero values */ 4095 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4096 4097 /* set matrix "i" values */ 4098 a->i[0] = 0; 4099 for (i=1; i<= M; i++) { 4100 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4101 a->ilen[i-1] = rowlengths[i-1]; 4102 } 4103 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4104 4105 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4106 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4107 PetscFunctionReturn(0); 4108 } 4109 4110 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4111 { 4112 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*b = (Mat_SeqAIJ*)B->data; 4113 PetscErrorCode ierr; 4114 #if defined(PETSC_USE_COMPLEX) 4115 PetscInt k; 4116 #endif 4117 4118 PetscFunctionBegin; 4119 /* If the matrix dimensions are not equal,or no of nonzeros */ 4120 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4121 *flg = PETSC_FALSE; 4122 PetscFunctionReturn(0); 4123 } 4124 4125 /* if the a->i are the same */ 4126 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4127 if (!*flg) PetscFunctionReturn(0); 4128 4129 /* if a->j are the same */ 4130 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4131 if (!*flg) PetscFunctionReturn(0); 4132 4133 /* if a->a are the same */ 4134 #if defined(PETSC_USE_COMPLEX) 4135 for (k=0; k<a->nz; k++) { 4136 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])) { 4137 *flg = PETSC_FALSE; 4138 PetscFunctionReturn(0); 4139 } 4140 } 4141 #else 4142 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4143 #endif 4144 PetscFunctionReturn(0); 4145 } 4146 4147 /*@ 4148 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4149 provided by the user. 4150 4151 Collective on MPI_Comm 4152 4153 Input Parameters: 4154 + comm - must be an MPI communicator of size 1 4155 . m - number of rows 4156 . n - number of columns 4157 . i - row indices 4158 . j - column indices 4159 - a - matrix values 4160 4161 Output Parameter: 4162 . mat - the matrix 4163 4164 Level: intermediate 4165 4166 Notes: 4167 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4168 once the matrix is destroyed and not before 4169 4170 You cannot set new nonzero locations into this matrix, that will generate an error. 4171 4172 The i and j indices are 0 based 4173 4174 The format which is used for the sparse matrix input, is equivalent to a 4175 row-major ordering.. i.e for the following matrix, the input data expected is 4176 as shown 4177 4178 $ 1 0 0 4179 $ 2 0 3 4180 $ 4 5 6 4181 $ 4182 $ i = {0,1,3,6} [size = nrow+1 = 3+1] 4183 $ j = {0,0,2,0,1,2} [size = 6]; values must be sorted for each row 4184 $ v = {1,2,3,4,5,6} [size = 6] 4185 4186 4187 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4188 4189 @*/ 4190 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat) 4191 { 4192 PetscErrorCode ierr; 4193 PetscInt ii; 4194 Mat_SeqAIJ *aij; 4195 #if defined(PETSC_USE_DEBUG) 4196 PetscInt jj; 4197 #endif 4198 4199 PetscFunctionBegin; 4200 if (m > 0 && i[0]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4201 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4202 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4203 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4204 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4205 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4206 aij = (Mat_SeqAIJ*)(*mat)->data; 4207 ierr = PetscMalloc2(m,&aij->imax,m,&aij->ilen);CHKERRQ(ierr); 4208 4209 aij->i = i; 4210 aij->j = j; 4211 aij->a = a; 4212 aij->singlemalloc = PETSC_FALSE; 4213 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4214 aij->free_a = PETSC_FALSE; 4215 aij->free_ij = PETSC_FALSE; 4216 4217 for (ii=0; ii<m; ii++) { 4218 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4219 #if defined(PETSC_USE_DEBUG) 4220 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]); 4221 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 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 not sorted",jj-i[ii],j[jj],ii); 4223 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); 4224 } 4225 #endif 4226 } 4227 #if defined(PETSC_USE_DEBUG) 4228 for (ii=0; ii<aij->i[m]; ii++) { 4229 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %D index = %D",ii,j[ii]); 4230 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]); 4231 } 4232 #endif 4233 4234 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4235 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4236 PetscFunctionReturn(0); 4237 } 4238 /*@C 4239 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4240 provided by the user. 4241 4242 Collective on MPI_Comm 4243 4244 Input Parameters: 4245 + comm - must be an MPI communicator of size 1 4246 . m - number of rows 4247 . n - number of columns 4248 . i - row indices 4249 . j - column indices 4250 . a - matrix values 4251 . nz - number of nonzeros 4252 - idx - 0 or 1 based 4253 4254 Output Parameter: 4255 . mat - the matrix 4256 4257 Level: intermediate 4258 4259 Notes: 4260 The i and j indices are 0 based 4261 4262 The format which is used for the sparse matrix input, is equivalent to a 4263 row-major ordering.. i.e for the following matrix, the input data expected is 4264 as shown: 4265 4266 1 0 0 4267 2 0 3 4268 4 5 6 4269 4270 i = {0,1,1,2,2,2} 4271 j = {0,0,2,0,1,2} 4272 v = {1,2,3,4,5,6} 4273 4274 4275 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4276 4277 @*/ 4278 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt i[],PetscInt j[],PetscScalar a[],Mat *mat,PetscInt nz,PetscBool idx) 4279 { 4280 PetscErrorCode ierr; 4281 PetscInt ii, *nnz, one = 1,row,col; 4282 4283 4284 PetscFunctionBegin; 4285 ierr = PetscCalloc1(m,&nnz);CHKERRQ(ierr); 4286 for (ii = 0; ii < nz; ii++) { 4287 nnz[i[ii] - !!idx] += 1; 4288 } 4289 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4290 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4291 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4292 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4293 for (ii = 0; ii < nz; ii++) { 4294 if (idx) { 4295 row = i[ii] - 1; 4296 col = j[ii] - 1; 4297 } else { 4298 row = i[ii]; 4299 col = j[ii]; 4300 } 4301 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4302 } 4303 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4304 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4305 ierr = PetscFree(nnz);CHKERRQ(ierr); 4306 PetscFunctionReturn(0); 4307 } 4308 4309 PetscErrorCode MatSeqAIJInvalidateDiagonal(Mat A) 4310 { 4311 Mat_SeqAIJ *a=(Mat_SeqAIJ*)A->data; 4312 PetscErrorCode ierr; 4313 4314 PetscFunctionBegin; 4315 a->idiagvalid = PETSC_FALSE; 4316 a->ibdiagvalid = PETSC_FALSE; 4317 4318 ierr = MatSeqAIJInvalidateDiagonal_Inode(A);CHKERRQ(ierr); 4319 PetscFunctionReturn(0); 4320 } 4321 4322 PetscErrorCode MatCreateMPIMatConcatenateSeqMat_SeqAIJ(MPI_Comm comm,Mat inmat,PetscInt n,MatReuse scall,Mat *outmat) 4323 { 4324 PetscErrorCode ierr; 4325 PetscMPIInt size; 4326 4327 PetscFunctionBegin; 4328 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4329 if (size == 1 && scall == MAT_REUSE_MATRIX) { 4330 ierr = MatCopy(inmat,*outmat,SAME_NONZERO_PATTERN);CHKERRQ(ierr); 4331 } else { 4332 ierr = MatCreateMPIMatConcatenateSeqMat_MPIAIJ(comm,inmat,n,scall,outmat);CHKERRQ(ierr); 4333 } 4334 PetscFunctionReturn(0); 4335 } 4336 4337 /* 4338 Permute A into C's *local* index space using rowemb,colemb. 4339 The embedding are supposed to be injections and the above implies that the range of rowemb is a subset 4340 of [0,m), colemb is in [0,n). 4341 If pattern == DIFFERENT_NONZERO_PATTERN, C is preallocated according to A. 4342 */ 4343 PetscErrorCode MatSetSeqMat_SeqAIJ(Mat C,IS rowemb,IS colemb,MatStructure pattern,Mat B) 4344 { 4345 /* If making this function public, change the error returned in this function away from _PLIB. */ 4346 PetscErrorCode ierr; 4347 Mat_SeqAIJ *Baij; 4348 PetscBool seqaij; 4349 PetscInt m,n,*nz,i,j,count; 4350 PetscScalar v; 4351 const PetscInt *rowindices,*colindices; 4352 4353 PetscFunctionBegin; 4354 if (!B) PetscFunctionReturn(0); 4355 /* Check to make sure the target matrix (and embeddings) are compatible with C and each other. */ 4356 ierr = PetscObjectTypeCompare((PetscObject)B,MATSEQAIJ,&seqaij);CHKERRQ(ierr); 4357 if (!seqaij) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is of wrong type"); 4358 if (rowemb) { 4359 ierr = ISGetLocalSize(rowemb,&m);CHKERRQ(ierr); 4360 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); 4361 } else { 4362 if (C->rmap->n != B->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is row-incompatible with the target matrix"); 4363 } 4364 if (colemb) { 4365 ierr = ISGetLocalSize(colemb,&n);CHKERRQ(ierr); 4366 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); 4367 } else { 4368 if (C->cmap->n != B->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Input matrix is col-incompatible with the target matrix"); 4369 } 4370 4371 Baij = (Mat_SeqAIJ*)(B->data); 4372 if (pattern == DIFFERENT_NONZERO_PATTERN) { 4373 ierr = PetscMalloc1(B->rmap->n,&nz);CHKERRQ(ierr); 4374 for (i=0; i<B->rmap->n; i++) { 4375 nz[i] = Baij->i[i+1] - Baij->i[i]; 4376 } 4377 ierr = MatSeqAIJSetPreallocation(C,0,nz);CHKERRQ(ierr); 4378 ierr = PetscFree(nz);CHKERRQ(ierr); 4379 } 4380 if (pattern == SUBSET_NONZERO_PATTERN) { 4381 ierr = MatZeroEntries(C);CHKERRQ(ierr); 4382 } 4383 count = 0; 4384 rowindices = NULL; 4385 colindices = NULL; 4386 if (rowemb) { 4387 ierr = ISGetIndices(rowemb,&rowindices);CHKERRQ(ierr); 4388 } 4389 if (colemb) { 4390 ierr = ISGetIndices(colemb,&colindices);CHKERRQ(ierr); 4391 } 4392 for (i=0; i<B->rmap->n; i++) { 4393 PetscInt row; 4394 row = i; 4395 if (rowindices) row = rowindices[i]; 4396 for (j=Baij->i[i]; j<Baij->i[i+1]; j++) { 4397 PetscInt col; 4398 col = Baij->j[count]; 4399 if (colindices) col = colindices[col]; 4400 v = Baij->a[count]; 4401 ierr = MatSetValues(C,1,&row,1,&col,&v,INSERT_VALUES);CHKERRQ(ierr); 4402 ++count; 4403 } 4404 } 4405 /* FIXME: set C's nonzerostate correctly. */ 4406 /* Assembly for C is necessary. */ 4407 C->preallocated = PETSC_TRUE; 4408 C->assembled = PETSC_TRUE; 4409 C->was_assembled = PETSC_FALSE; 4410 PetscFunctionReturn(0); 4411 } 4412 4413 4414 /* 4415 Special version for direct calls from Fortran 4416 */ 4417 #include <petsc/private/fortranimpl.h> 4418 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4419 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4420 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4421 #define matsetvaluesseqaij_ matsetvaluesseqaij 4422 #endif 4423 4424 /* Change these macros so can be used in void function */ 4425 #undef CHKERRQ 4426 #define CHKERRQ(ierr) CHKERRABORT(PetscObjectComm((PetscObject)A),ierr) 4427 #undef SETERRQ2 4428 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4429 #undef SETERRQ3 4430 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4431 4432 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) 4433 { 4434 Mat A = *AA; 4435 PetscInt m = *mm, n = *nn; 4436 InsertMode is = *isis; 4437 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4438 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4439 PetscInt *imax,*ai,*ailen; 4440 PetscErrorCode ierr; 4441 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4442 MatScalar *ap,value,*aa; 4443 PetscBool ignorezeroentries = a->ignorezeroentries; 4444 PetscBool roworiented = a->roworiented; 4445 4446 PetscFunctionBegin; 4447 MatCheckPreallocated(A,1); 4448 imax = a->imax; 4449 ai = a->i; 4450 ailen = a->ilen; 4451 aj = a->j; 4452 aa = a->a; 4453 4454 for (k=0; k<m; k++) { /* loop over added rows */ 4455 row = im[k]; 4456 if (row < 0) continue; 4457 #if defined(PETSC_USE_DEBUG) 4458 if (row >= A->rmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4459 #endif 4460 rp = aj + ai[row]; ap = aa + ai[row]; 4461 rmax = imax[row]; nrow = ailen[row]; 4462 low = 0; 4463 high = nrow; 4464 for (l=0; l<n; l++) { /* loop over added columns */ 4465 if (in[l] < 0) continue; 4466 #if defined(PETSC_USE_DEBUG) 4467 if (in[l] >= A->cmap->n) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4468 #endif 4469 col = in[l]; 4470 if (roworiented) value = v[l + k*n]; 4471 else value = v[k + l*m]; 4472 4473 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4474 4475 if (col <= lastcol) low = 0; 4476 else high = nrow; 4477 lastcol = col; 4478 while (high-low > 5) { 4479 t = (low+high)/2; 4480 if (rp[t] > col) high = t; 4481 else low = t; 4482 } 4483 for (i=low; i<high; i++) { 4484 if (rp[i] > col) break; 4485 if (rp[i] == col) { 4486 if (is == ADD_VALUES) ap[i] += value; 4487 else ap[i] = value; 4488 goto noinsert; 4489 } 4490 } 4491 if (value == 0.0 && ignorezeroentries) goto noinsert; 4492 if (nonew == 1) goto noinsert; 4493 if (nonew == -1) SETERRABORT(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4494 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4495 N = nrow++ - 1; a->nz++; high++; 4496 /* shift up all the later entries in this row */ 4497 for (ii=N; ii>=i; ii--) { 4498 rp[ii+1] = rp[ii]; 4499 ap[ii+1] = ap[ii]; 4500 } 4501 rp[i] = col; 4502 ap[i] = value; 4503 A->nonzerostate++; 4504 noinsert:; 4505 low = i + 1; 4506 } 4507 ailen[row] = nrow; 4508 } 4509 PetscFunctionReturnVoid(); 4510 } 4511 4512