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