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