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