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