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