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