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 sum += PetscRealPart(PetscConj(*v)*(*v)); v++; 1809 } 1810 *nrm = PetscSqrtReal(sum); 1811 } else if (type == NORM_1) { 1812 PetscReal *tmp; 1813 PetscInt *jj = a->j; 1814 ierr = PetscMalloc((A->cmap->n+1)*sizeof(PetscReal),&tmp);CHKERRQ(ierr); 1815 ierr = PetscMemzero(tmp,A->cmap->n*sizeof(PetscReal));CHKERRQ(ierr); 1816 *nrm = 0.0; 1817 for (j=0; j<a->nz; j++) { 1818 tmp[*jj++] += PetscAbsScalar(*v); v++; 1819 } 1820 for (j=0; j<A->cmap->n; j++) { 1821 if (tmp[j] > *nrm) *nrm = tmp[j]; 1822 } 1823 ierr = PetscFree(tmp);CHKERRQ(ierr); 1824 } else if (type == NORM_INFINITY) { 1825 *nrm = 0.0; 1826 for (j=0; j<A->rmap->n; j++) { 1827 v = a->a + a->i[j]; 1828 sum = 0.0; 1829 for (i=0; i<a->i[j+1]-a->i[j]; i++) { 1830 sum += PetscAbsScalar(*v); v++; 1831 } 1832 if (sum > *nrm) *nrm = sum; 1833 } 1834 } else { 1835 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No support for two norm"); 1836 } 1837 PetscFunctionReturn(0); 1838 } 1839 1840 /* Merged from MatGetSymbolicTranspose_SeqAIJ() - replace MatGetSymbolicTranspose_SeqAIJ()? */ 1841 #undef __FUNCT__ 1842 #define __FUNCT__ "MatTransposeSymbolic_SeqAIJ" 1843 PetscErrorCode MatTransposeSymbolic_SeqAIJ(Mat A,Mat *B) 1844 { 1845 PetscErrorCode ierr; 1846 PetscInt i,j,anzj; 1847 Mat_SeqAIJ *a=(Mat_SeqAIJ *)A->data,*b; 1848 PetscInt an=A->cmap->N,am=A->rmap->N; 1849 PetscInt *ati,*atj,*atfill,*ai=a->i,*aj=a->j; 1850 1851 PetscFunctionBegin; 1852 /* Allocate space for symbolic transpose info and work array */ 1853 ierr = PetscMalloc((an+1)*sizeof(PetscInt),&ati);CHKERRQ(ierr); 1854 ierr = PetscMalloc(ai[am]*sizeof(PetscInt),&atj);CHKERRQ(ierr); 1855 ierr = PetscMalloc(an*sizeof(PetscInt),&atfill);CHKERRQ(ierr); 1856 ierr = PetscMemzero(ati,(an+1)*sizeof(PetscInt));CHKERRQ(ierr); 1857 1858 /* Walk through aj and count ## of non-zeros in each row of A^T. */ 1859 /* Note: offset by 1 for fast conversion into csr format. */ 1860 for (i=0;i<ai[am];i++) { 1861 ati[aj[i]+1] += 1; 1862 } 1863 /* Form ati for csr format of A^T. */ 1864 for (i=0;i<an;i++) { 1865 ati[i+1] += ati[i]; 1866 } 1867 1868 /* Copy ati into atfill so we have locations of the next free space in atj */ 1869 ierr = PetscMemcpy(atfill,ati,an*sizeof(PetscInt));CHKERRQ(ierr); 1870 1871 /* Walk through A row-wise and mark nonzero entries of A^T. */ 1872 for (i=0;i<am;i++) { 1873 anzj = ai[i+1] - ai[i]; 1874 for (j=0;j<anzj;j++) { 1875 atj[atfill[*aj]] = i; 1876 atfill[*aj++] += 1; 1877 } 1878 } 1879 1880 /* Clean up temporary space and complete requests. */ 1881 ierr = PetscFree(atfill);CHKERRQ(ierr); 1882 ierr = MatCreateSeqAIJWithArrays(((PetscObject)A)->comm,an,am,ati,atj,PETSC_NULL,B);CHKERRQ(ierr); 1883 (*B)->rmap->bs = A->cmap->bs; 1884 (*B)->cmap->bs = A->rmap->bs; 1885 1886 b = (Mat_SeqAIJ *)((*B)->data); 1887 b->free_a = PETSC_FALSE; 1888 b->free_ij = PETSC_TRUE; 1889 b->nonew = 0; 1890 PetscFunctionReturn(0); 1891 } 1892 1893 #undef __FUNCT__ 1894 #define __FUNCT__ "MatTranspose_SeqAIJ" 1895 PetscErrorCode MatTranspose_SeqAIJ(Mat A,MatReuse reuse,Mat *B) 1896 { 1897 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 1898 Mat C; 1899 PetscErrorCode ierr; 1900 PetscInt i,*aj = a->j,*ai = a->i,m = A->rmap->n,len,*col; 1901 MatScalar *array = a->a; 1902 1903 PetscFunctionBegin; 1904 if (reuse == MAT_REUSE_MATRIX && A == *B && m != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Square matrix only for in-place"); 1905 1906 if (reuse == MAT_INITIAL_MATRIX || *B == A) { 1907 ierr = PetscMalloc((1+A->cmap->n)*sizeof(PetscInt),&col);CHKERRQ(ierr); 1908 ierr = PetscMemzero(col,(1+A->cmap->n)*sizeof(PetscInt));CHKERRQ(ierr); 1909 1910 for (i=0; i<ai[m]; i++) col[aj[i]] += 1; 1911 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 1912 ierr = MatSetSizes(C,A->cmap->n,m,A->cmap->n,m);CHKERRQ(ierr); 1913 ierr = MatSetBlockSizes(C,A->cmap->bs,A->rmap->bs);CHKERRQ(ierr); 1914 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 1915 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,col);CHKERRQ(ierr); 1916 ierr = PetscFree(col);CHKERRQ(ierr); 1917 } else { 1918 C = *B; 1919 } 1920 1921 for (i=0; i<m; i++) { 1922 len = ai[i+1]-ai[i]; 1923 ierr = MatSetValues_SeqAIJ(C,len,aj,1,&i,array,INSERT_VALUES);CHKERRQ(ierr); 1924 array += len; 1925 aj += len; 1926 } 1927 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1928 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 1929 1930 if (reuse == MAT_INITIAL_MATRIX || *B != A) { 1931 *B = C; 1932 } else { 1933 ierr = MatHeaderMerge(A,C);CHKERRQ(ierr); 1934 } 1935 PetscFunctionReturn(0); 1936 } 1937 1938 EXTERN_C_BEGIN 1939 #undef __FUNCT__ 1940 #define __FUNCT__ "MatIsTranspose_SeqAIJ" 1941 PetscErrorCode MatIsTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 1942 { 1943 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1944 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1945 MatScalar *va,*vb; 1946 PetscErrorCode ierr; 1947 PetscInt ma,na,mb,nb, i; 1948 1949 PetscFunctionBegin; 1950 bij = (Mat_SeqAIJ *) B->data; 1951 1952 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 1953 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 1954 if (ma!=nb || na!=mb){ 1955 *f = PETSC_FALSE; 1956 PetscFunctionReturn(0); 1957 } 1958 aii = aij->i; bii = bij->i; 1959 adx = aij->j; bdx = bij->j; 1960 va = aij->a; vb = bij->a; 1961 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 1962 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 1963 for (i=0; i<ma; i++) aptr[i] = aii[i]; 1964 for (i=0; i<mb; i++) bptr[i] = bii[i]; 1965 1966 *f = PETSC_TRUE; 1967 for (i=0; i<ma; i++) { 1968 while (aptr[i]<aii[i+1]) { 1969 PetscInt idc,idr; 1970 PetscScalar vc,vr; 1971 /* column/row index/value */ 1972 idc = adx[aptr[i]]; 1973 idr = bdx[bptr[idc]]; 1974 vc = va[aptr[i]]; 1975 vr = vb[bptr[idc]]; 1976 if (i!=idr || PetscAbsScalar(vc-vr) > tol) { 1977 *f = PETSC_FALSE; 1978 goto done; 1979 } else { 1980 aptr[i]++; 1981 if (B || i!=idc) bptr[idc]++; 1982 } 1983 } 1984 } 1985 done: 1986 ierr = PetscFree(aptr);CHKERRQ(ierr); 1987 ierr = PetscFree(bptr);CHKERRQ(ierr); 1988 PetscFunctionReturn(0); 1989 } 1990 EXTERN_C_END 1991 1992 EXTERN_C_BEGIN 1993 #undef __FUNCT__ 1994 #define __FUNCT__ "MatIsHermitianTranspose_SeqAIJ" 1995 PetscErrorCode MatIsHermitianTranspose_SeqAIJ(Mat A,Mat B,PetscReal tol,PetscBool *f) 1996 { 1997 Mat_SeqAIJ *aij = (Mat_SeqAIJ *) A->data,*bij = (Mat_SeqAIJ*) A->data; 1998 PetscInt *adx,*bdx,*aii,*bii,*aptr,*bptr; 1999 MatScalar *va,*vb; 2000 PetscErrorCode ierr; 2001 PetscInt ma,na,mb,nb, i; 2002 2003 PetscFunctionBegin; 2004 bij = (Mat_SeqAIJ *) B->data; 2005 2006 ierr = MatGetSize(A,&ma,&na);CHKERRQ(ierr); 2007 ierr = MatGetSize(B,&mb,&nb);CHKERRQ(ierr); 2008 if (ma!=nb || na!=mb){ 2009 *f = PETSC_FALSE; 2010 PetscFunctionReturn(0); 2011 } 2012 aii = aij->i; bii = bij->i; 2013 adx = aij->j; bdx = bij->j; 2014 va = aij->a; vb = bij->a; 2015 ierr = PetscMalloc(ma*sizeof(PetscInt),&aptr);CHKERRQ(ierr); 2016 ierr = PetscMalloc(mb*sizeof(PetscInt),&bptr);CHKERRQ(ierr); 2017 for (i=0; i<ma; i++) aptr[i] = aii[i]; 2018 for (i=0; i<mb; i++) bptr[i] = bii[i]; 2019 2020 *f = PETSC_TRUE; 2021 for (i=0; i<ma; i++) { 2022 while (aptr[i]<aii[i+1]) { 2023 PetscInt idc,idr; 2024 PetscScalar vc,vr; 2025 /* column/row index/value */ 2026 idc = adx[aptr[i]]; 2027 idr = bdx[bptr[idc]]; 2028 vc = va[aptr[i]]; 2029 vr = vb[bptr[idc]]; 2030 if (i!=idr || PetscAbsScalar(vc-PetscConj(vr)) > tol) { 2031 *f = PETSC_FALSE; 2032 goto done; 2033 } else { 2034 aptr[i]++; 2035 if (B || i!=idc) bptr[idc]++; 2036 } 2037 } 2038 } 2039 done: 2040 ierr = PetscFree(aptr);CHKERRQ(ierr); 2041 ierr = PetscFree(bptr);CHKERRQ(ierr); 2042 PetscFunctionReturn(0); 2043 } 2044 EXTERN_C_END 2045 2046 #undef __FUNCT__ 2047 #define __FUNCT__ "MatIsSymmetric_SeqAIJ" 2048 PetscErrorCode MatIsSymmetric_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2049 { 2050 PetscErrorCode ierr; 2051 PetscFunctionBegin; 2052 ierr = MatIsTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2053 PetscFunctionReturn(0); 2054 } 2055 2056 #undef __FUNCT__ 2057 #define __FUNCT__ "MatIsHermitian_SeqAIJ" 2058 PetscErrorCode MatIsHermitian_SeqAIJ(Mat A,PetscReal tol,PetscBool *f) 2059 { 2060 PetscErrorCode ierr; 2061 PetscFunctionBegin; 2062 ierr = MatIsHermitianTranspose_SeqAIJ(A,A,tol,f);CHKERRQ(ierr); 2063 PetscFunctionReturn(0); 2064 } 2065 2066 #undef __FUNCT__ 2067 #define __FUNCT__ "MatDiagonalScale_SeqAIJ" 2068 PetscErrorCode MatDiagonalScale_SeqAIJ(Mat A,Vec ll,Vec rr) 2069 { 2070 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2071 PetscScalar *l,*r,x; 2072 MatScalar *v; 2073 PetscErrorCode ierr; 2074 PetscInt i,j,m = A->rmap->n,n = A->cmap->n,M,nz = a->nz,*jj; 2075 2076 PetscFunctionBegin; 2077 if (ll) { 2078 /* The local size is used so that VecMPI can be passed to this routine 2079 by MatDiagonalScale_MPIAIJ */ 2080 ierr = VecGetLocalSize(ll,&m);CHKERRQ(ierr); 2081 if (m != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Left scaling vector wrong length"); 2082 ierr = VecGetArray(ll,&l);CHKERRQ(ierr); 2083 v = a->a; 2084 for (i=0; i<m; i++) { 2085 x = l[i]; 2086 M = a->i[i+1] - a->i[i]; 2087 for (j=0; j<M; j++) { (*v++) *= x;} 2088 } 2089 ierr = VecRestoreArray(ll,&l);CHKERRQ(ierr); 2090 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2091 } 2092 if (rr) { 2093 ierr = VecGetLocalSize(rr,&n);CHKERRQ(ierr); 2094 if (n != A->cmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Right scaling vector wrong length"); 2095 ierr = VecGetArray(rr,&r);CHKERRQ(ierr); 2096 v = a->a; jj = a->j; 2097 for (i=0; i<nz; i++) { 2098 (*v++) *= r[*jj++]; 2099 } 2100 ierr = VecRestoreArray(rr,&r);CHKERRQ(ierr); 2101 ierr = PetscLogFlops(nz);CHKERRQ(ierr); 2102 } 2103 a->idiagvalid = PETSC_FALSE; 2104 a->ibdiagvalid = PETSC_FALSE; 2105 PetscFunctionReturn(0); 2106 } 2107 2108 #undef __FUNCT__ 2109 #define __FUNCT__ "MatGetSubMatrix_SeqAIJ" 2110 PetscErrorCode MatGetSubMatrix_SeqAIJ(Mat A,IS isrow,IS iscol,PetscInt csize,MatReuse scall,Mat *B) 2111 { 2112 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data,*c; 2113 PetscErrorCode ierr; 2114 PetscInt *smap,i,k,kstart,kend,oldcols = A->cmap->n,*lens; 2115 PetscInt row,mat_i,*mat_j,tcol,first,step,*mat_ilen,sum,lensi; 2116 const PetscInt *irow,*icol; 2117 PetscInt nrows,ncols; 2118 PetscInt *starts,*j_new,*i_new,*aj = a->j,*ai = a->i,ii,*ailen = a->ilen; 2119 MatScalar *a_new,*mat_a; 2120 Mat C; 2121 PetscBool stride,sorted; 2122 2123 PetscFunctionBegin; 2124 ierr = ISSorted(isrow,&sorted);CHKERRQ(ierr); 2125 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"ISrow is not sorted"); 2126 ierr = ISSorted(iscol,&sorted);CHKERRQ(ierr); 2127 if (!sorted) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"IScol is not sorted"); 2128 2129 ierr = ISGetIndices(isrow,&irow);CHKERRQ(ierr); 2130 ierr = ISGetLocalSize(isrow,&nrows);CHKERRQ(ierr); 2131 ierr = ISGetLocalSize(iscol,&ncols);CHKERRQ(ierr); 2132 2133 ierr = ISStrideGetInfo(iscol,&first,&step);CHKERRQ(ierr); 2134 ierr = PetscObjectTypeCompare((PetscObject)iscol,ISSTRIDE,&stride);CHKERRQ(ierr); 2135 if (stride && step == 1) { 2136 /* special case of contiguous rows */ 2137 ierr = PetscMalloc2(nrows,PetscInt,&lens,nrows,PetscInt,&starts);CHKERRQ(ierr); 2138 /* loop over new rows determining lens and starting points */ 2139 for (i=0; i<nrows; i++) { 2140 kstart = ai[irow[i]]; 2141 kend = kstart + ailen[irow[i]]; 2142 for (k=kstart; k<kend; k++) { 2143 if (aj[k] >= first) { 2144 starts[i] = k; 2145 break; 2146 } 2147 } 2148 sum = 0; 2149 while (k < kend) { 2150 if (aj[k++] >= first+ncols) break; 2151 sum++; 2152 } 2153 lens[i] = sum; 2154 } 2155 /* create submatrix */ 2156 if (scall == MAT_REUSE_MATRIX) { 2157 PetscInt n_cols,n_rows; 2158 ierr = MatGetSize(*B,&n_rows,&n_cols);CHKERRQ(ierr); 2159 if (n_rows != nrows || n_cols != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Reused submatrix wrong size"); 2160 ierr = MatZeroEntries(*B);CHKERRQ(ierr); 2161 C = *B; 2162 } else { 2163 PetscInt rbs,cbs; 2164 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2165 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2166 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2167 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2168 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2169 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2170 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2171 } 2172 c = (Mat_SeqAIJ*)C->data; 2173 2174 /* loop over rows inserting into submatrix */ 2175 a_new = c->a; 2176 j_new = c->j; 2177 i_new = c->i; 2178 2179 for (i=0; i<nrows; i++) { 2180 ii = starts[i]; 2181 lensi = lens[i]; 2182 for (k=0; k<lensi; k++) { 2183 *j_new++ = aj[ii+k] - first; 2184 } 2185 ierr = PetscMemcpy(a_new,a->a + starts[i],lensi*sizeof(PetscScalar));CHKERRQ(ierr); 2186 a_new += lensi; 2187 i_new[i+1] = i_new[i] + lensi; 2188 c->ilen[i] = lensi; 2189 } 2190 ierr = PetscFree2(lens,starts);CHKERRQ(ierr); 2191 } else { 2192 ierr = ISGetIndices(iscol,&icol);CHKERRQ(ierr); 2193 ierr = PetscMalloc(oldcols*sizeof(PetscInt),&smap);CHKERRQ(ierr); 2194 ierr = PetscMemzero(smap,oldcols*sizeof(PetscInt));CHKERRQ(ierr); 2195 ierr = PetscMalloc((1+nrows)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2196 for (i=0; i<ncols; i++) { 2197 #if defined(PETSC_USE_DEBUG) 2198 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); 2199 #endif 2200 smap[icol[i]] = i+1; 2201 } 2202 2203 /* determine lens of each row */ 2204 for (i=0; i<nrows; i++) { 2205 kstart = ai[irow[i]]; 2206 kend = kstart + a->ilen[irow[i]]; 2207 lens[i] = 0; 2208 for (k=kstart; k<kend; k++) { 2209 if (smap[aj[k]]) { 2210 lens[i]++; 2211 } 2212 } 2213 } 2214 /* Create and fill new matrix */ 2215 if (scall == MAT_REUSE_MATRIX) { 2216 PetscBool equal; 2217 2218 c = (Mat_SeqAIJ *)((*B)->data); 2219 if ((*B)->rmap->n != nrows || (*B)->cmap->n != ncols) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong size"); 2220 ierr = PetscMemcmp(c->ilen,lens,(*B)->rmap->n*sizeof(PetscInt),&equal);CHKERRQ(ierr); 2221 if (!equal) { 2222 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Cannot reuse matrix. wrong no of nonzeros"); 2223 } 2224 ierr = PetscMemzero(c->ilen,(*B)->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 2225 C = *B; 2226 } else { 2227 PetscInt rbs,cbs; 2228 ierr = MatCreate(((PetscObject)A)->comm,&C);CHKERRQ(ierr); 2229 ierr = MatSetSizes(C,nrows,ncols,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 2230 ierr = ISGetBlockSize(isrow,&rbs);CHKERRQ(ierr); 2231 ierr = ISGetBlockSize(iscol,&cbs);CHKERRQ(ierr); 2232 ierr = MatSetBlockSizes(C,rbs,cbs);CHKERRQ(ierr); 2233 ierr = MatSetType(C,((PetscObject)A)->type_name);CHKERRQ(ierr); 2234 ierr = MatSeqAIJSetPreallocation_SeqAIJ(C,0,lens);CHKERRQ(ierr); 2235 } 2236 c = (Mat_SeqAIJ *)(C->data); 2237 for (i=0; i<nrows; i++) { 2238 row = irow[i]; 2239 kstart = ai[row]; 2240 kend = kstart + a->ilen[row]; 2241 mat_i = c->i[i]; 2242 mat_j = c->j + mat_i; 2243 mat_a = c->a + mat_i; 2244 mat_ilen = c->ilen + i; 2245 for (k=kstart; k<kend; k++) { 2246 if ((tcol=smap[a->j[k]])) { 2247 *mat_j++ = tcol - 1; 2248 *mat_a++ = a->a[k]; 2249 (*mat_ilen)++; 2250 2251 } 2252 } 2253 } 2254 /* Free work space */ 2255 ierr = ISRestoreIndices(iscol,&icol);CHKERRQ(ierr); 2256 ierr = PetscFree(smap);CHKERRQ(ierr); 2257 ierr = PetscFree(lens);CHKERRQ(ierr); 2258 } 2259 ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2260 ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2261 2262 ierr = ISRestoreIndices(isrow,&irow);CHKERRQ(ierr); 2263 *B = C; 2264 PetscFunctionReturn(0); 2265 } 2266 2267 #undef __FUNCT__ 2268 #define __FUNCT__ "MatGetMultiProcBlock_SeqAIJ" 2269 PetscErrorCode MatGetMultiProcBlock_SeqAIJ(Mat mat,MPI_Comm subComm,MatReuse scall,Mat* subMat) 2270 { 2271 PetscErrorCode ierr; 2272 Mat B; 2273 2274 PetscFunctionBegin; 2275 ierr = MatCreate(subComm,&B);CHKERRQ(ierr); 2276 ierr = MatSetSizes(B,mat->rmap->n,mat->cmap->n,mat->rmap->n,mat->cmap->n);CHKERRQ(ierr); 2277 ierr = MatSetBlockSizes(B,mat->rmap->bs,mat->cmap->bs);CHKERRQ(ierr); 2278 ierr = MatSetType(B,MATSEQAIJ);CHKERRQ(ierr); 2279 ierr = MatDuplicateNoCreate_SeqAIJ(B,mat,MAT_COPY_VALUES,PETSC_TRUE);CHKERRQ(ierr); 2280 *subMat = B; 2281 PetscFunctionReturn(0); 2282 } 2283 2284 #undef __FUNCT__ 2285 #define __FUNCT__ "MatILUFactor_SeqAIJ" 2286 PetscErrorCode MatILUFactor_SeqAIJ(Mat inA,IS row,IS col,const MatFactorInfo *info) 2287 { 2288 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2289 PetscErrorCode ierr; 2290 Mat outA; 2291 PetscBool row_identity,col_identity; 2292 2293 PetscFunctionBegin; 2294 if (info->levels != 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Only levels=0 supported for in-place ilu"); 2295 2296 ierr = ISIdentity(row,&row_identity);CHKERRQ(ierr); 2297 ierr = ISIdentity(col,&col_identity);CHKERRQ(ierr); 2298 2299 outA = inA; 2300 outA->factortype = MAT_FACTOR_LU; 2301 ierr = PetscObjectReference((PetscObject)row);CHKERRQ(ierr); 2302 ierr = ISDestroy(&a->row);CHKERRQ(ierr); 2303 a->row = row; 2304 ierr = PetscObjectReference((PetscObject)col);CHKERRQ(ierr); 2305 ierr = ISDestroy(&a->col);CHKERRQ(ierr); 2306 a->col = col; 2307 2308 /* Create the inverse permutation so that it can be used in MatLUFactorNumeric() */ 2309 ierr = ISDestroy(&a->icol);CHKERRQ(ierr); 2310 ierr = ISInvertPermutation(col,PETSC_DECIDE,&a->icol);CHKERRQ(ierr); 2311 ierr = PetscLogObjectParent(inA,a->icol);CHKERRQ(ierr); 2312 2313 if (!a->solve_work) { /* this matrix may have been factored before */ 2314 ierr = PetscMalloc((inA->rmap->n+1)*sizeof(PetscScalar),&a->solve_work);CHKERRQ(ierr); 2315 ierr = PetscLogObjectMemory(inA, (inA->rmap->n+1)*sizeof(PetscScalar));CHKERRQ(ierr); 2316 } 2317 2318 ierr = MatMarkDiagonal_SeqAIJ(inA);CHKERRQ(ierr); 2319 if (row_identity && col_identity) { 2320 ierr = MatLUFactorNumeric_SeqAIJ_inplace(outA,inA,info);CHKERRQ(ierr); 2321 } else { 2322 ierr = MatLUFactorNumeric_SeqAIJ_InplaceWithPerm(outA,inA,info);CHKERRQ(ierr); 2323 } 2324 PetscFunctionReturn(0); 2325 } 2326 2327 #undef __FUNCT__ 2328 #define __FUNCT__ "MatScale_SeqAIJ" 2329 PetscErrorCode MatScale_SeqAIJ(Mat inA,PetscScalar alpha) 2330 { 2331 Mat_SeqAIJ *a = (Mat_SeqAIJ*)inA->data; 2332 PetscScalar oalpha = alpha; 2333 PetscErrorCode ierr; 2334 PetscBLASInt one = 1,bnz = PetscBLASIntCast(a->nz); 2335 2336 PetscFunctionBegin; 2337 BLASscal_(&bnz,&oalpha,a->a,&one); 2338 ierr = PetscLogFlops(a->nz);CHKERRQ(ierr); 2339 a->idiagvalid = PETSC_FALSE; 2340 a->ibdiagvalid = PETSC_FALSE; 2341 PetscFunctionReturn(0); 2342 } 2343 2344 #undef __FUNCT__ 2345 #define __FUNCT__ "MatGetSubMatrices_SeqAIJ" 2346 PetscErrorCode MatGetSubMatrices_SeqAIJ(Mat A,PetscInt n,const IS irow[],const IS icol[],MatReuse scall,Mat *B[]) 2347 { 2348 PetscErrorCode ierr; 2349 PetscInt i; 2350 2351 PetscFunctionBegin; 2352 if (scall == MAT_INITIAL_MATRIX) { 2353 ierr = PetscMalloc((n+1)*sizeof(Mat),B);CHKERRQ(ierr); 2354 } 2355 2356 for (i=0; i<n; i++) { 2357 ierr = MatGetSubMatrix_SeqAIJ(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(ierr); 2358 } 2359 PetscFunctionReturn(0); 2360 } 2361 2362 #undef __FUNCT__ 2363 #define __FUNCT__ "MatIncreaseOverlap_SeqAIJ" 2364 PetscErrorCode MatIncreaseOverlap_SeqAIJ(Mat A,PetscInt is_max,IS is[],PetscInt ov) 2365 { 2366 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2367 PetscErrorCode ierr; 2368 PetscInt row,i,j,k,l,m,n,*nidx,isz,val; 2369 const PetscInt *idx; 2370 PetscInt start,end,*ai,*aj; 2371 PetscBT table; 2372 2373 PetscFunctionBegin; 2374 m = A->rmap->n; 2375 ai = a->i; 2376 aj = a->j; 2377 2378 if (ov < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"illegal negative overlap value used"); 2379 2380 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&nidx);CHKERRQ(ierr); 2381 ierr = PetscBTCreate(m,&table);CHKERRQ(ierr); 2382 2383 for (i=0; i<is_max; i++) { 2384 /* Initialize the two local arrays */ 2385 isz = 0; 2386 ierr = PetscBTMemzero(m,table);CHKERRQ(ierr); 2387 2388 /* Extract the indices, assume there can be duplicate entries */ 2389 ierr = ISGetIndices(is[i],&idx);CHKERRQ(ierr); 2390 ierr = ISGetLocalSize(is[i],&n);CHKERRQ(ierr); 2391 2392 /* Enter these into the temp arrays. I.e., mark table[row], enter row into new index */ 2393 for (j=0; j<n ; ++j){ 2394 if (!PetscBTLookupSet(table,idx[j])) { nidx[isz++] = idx[j];} 2395 } 2396 ierr = ISRestoreIndices(is[i],&idx);CHKERRQ(ierr); 2397 ierr = ISDestroy(&is[i]);CHKERRQ(ierr); 2398 2399 k = 0; 2400 for (j=0; j<ov; j++){ /* for each overlap */ 2401 n = isz; 2402 for (; k<n ; k++){ /* do only those rows in nidx[k], which are not done yet */ 2403 row = nidx[k]; 2404 start = ai[row]; 2405 end = ai[row+1]; 2406 for (l = start; l<end ; l++){ 2407 val = aj[l] ; 2408 if (!PetscBTLookupSet(table,val)) {nidx[isz++] = val;} 2409 } 2410 } 2411 } 2412 ierr = ISCreateGeneral(PETSC_COMM_SELF,isz,nidx,PETSC_COPY_VALUES,(is+i));CHKERRQ(ierr); 2413 } 2414 ierr = PetscBTDestroy(&table);CHKERRQ(ierr); 2415 ierr = PetscFree(nidx);CHKERRQ(ierr); 2416 PetscFunctionReturn(0); 2417 } 2418 2419 /* -------------------------------------------------------------- */ 2420 #undef __FUNCT__ 2421 #define __FUNCT__ "MatPermute_SeqAIJ" 2422 PetscErrorCode MatPermute_SeqAIJ(Mat A,IS rowp,IS colp,Mat *B) 2423 { 2424 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2425 PetscErrorCode ierr; 2426 PetscInt i,nz = 0,m = A->rmap->n,n = A->cmap->n; 2427 const PetscInt *row,*col; 2428 PetscInt *cnew,j,*lens; 2429 IS icolp,irowp; 2430 PetscInt *cwork = PETSC_NULL; 2431 PetscScalar *vwork = PETSC_NULL; 2432 2433 PetscFunctionBegin; 2434 ierr = ISInvertPermutation(rowp,PETSC_DECIDE,&irowp);CHKERRQ(ierr); 2435 ierr = ISGetIndices(irowp,&row);CHKERRQ(ierr); 2436 ierr = ISInvertPermutation(colp,PETSC_DECIDE,&icolp);CHKERRQ(ierr); 2437 ierr = ISGetIndices(icolp,&col);CHKERRQ(ierr); 2438 2439 /* determine lengths of permuted rows */ 2440 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&lens);CHKERRQ(ierr); 2441 for (i=0; i<m; i++) { 2442 lens[row[i]] = a->i[i+1] - a->i[i]; 2443 } 2444 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 2445 ierr = MatSetSizes(*B,m,n,m,n);CHKERRQ(ierr); 2446 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 2447 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 2448 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*B,0,lens);CHKERRQ(ierr); 2449 ierr = PetscFree(lens);CHKERRQ(ierr); 2450 2451 ierr = PetscMalloc(n*sizeof(PetscInt),&cnew);CHKERRQ(ierr); 2452 for (i=0; i<m; i++) { 2453 ierr = MatGetRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2454 for (j=0; j<nz; j++) { cnew[j] = col[cwork[j]];} 2455 ierr = MatSetValues_SeqAIJ(*B,1,&row[i],nz,cnew,vwork,INSERT_VALUES);CHKERRQ(ierr); 2456 ierr = MatRestoreRow_SeqAIJ(A,i,&nz,&cwork,&vwork);CHKERRQ(ierr); 2457 } 2458 ierr = PetscFree(cnew);CHKERRQ(ierr); 2459 (*B)->assembled = PETSC_FALSE; 2460 ierr = MatAssemblyBegin(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2461 ierr = MatAssemblyEnd(*B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2462 ierr = ISRestoreIndices(irowp,&row);CHKERRQ(ierr); 2463 ierr = ISRestoreIndices(icolp,&col);CHKERRQ(ierr); 2464 ierr = ISDestroy(&irowp);CHKERRQ(ierr); 2465 ierr = ISDestroy(&icolp);CHKERRQ(ierr); 2466 PetscFunctionReturn(0); 2467 } 2468 2469 #undef __FUNCT__ 2470 #define __FUNCT__ "MatCopy_SeqAIJ" 2471 PetscErrorCode MatCopy_SeqAIJ(Mat A,Mat B,MatStructure str) 2472 { 2473 PetscErrorCode ierr; 2474 2475 PetscFunctionBegin; 2476 /* If the two matrices have the same copy implementation, use fast copy. */ 2477 if (str == SAME_NONZERO_PATTERN && (A->ops->copy == B->ops->copy)) { 2478 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2479 Mat_SeqAIJ *b = (Mat_SeqAIJ*)B->data; 2480 2481 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"); 2482 ierr = PetscMemcpy(b->a,a->a,(a->i[A->rmap->n])*sizeof(PetscScalar));CHKERRQ(ierr); 2483 } else { 2484 ierr = MatCopy_Basic(A,B,str);CHKERRQ(ierr); 2485 } 2486 PetscFunctionReturn(0); 2487 } 2488 2489 #undef __FUNCT__ 2490 #define __FUNCT__ "MatSetUp_SeqAIJ" 2491 PetscErrorCode MatSetUp_SeqAIJ(Mat A) 2492 { 2493 PetscErrorCode ierr; 2494 2495 PetscFunctionBegin; 2496 ierr = MatSeqAIJSetPreallocation_SeqAIJ(A,PETSC_DEFAULT,0);CHKERRQ(ierr); 2497 PetscFunctionReturn(0); 2498 } 2499 2500 #undef __FUNCT__ 2501 #define __FUNCT__ "MatSeqAIJGetArray_SeqAIJ" 2502 PetscErrorCode MatSeqAIJGetArray_SeqAIJ(Mat A,PetscScalar *array[]) 2503 { 2504 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2505 PetscFunctionBegin; 2506 *array = a->a; 2507 PetscFunctionReturn(0); 2508 } 2509 2510 #undef __FUNCT__ 2511 #define __FUNCT__ "MatSeqAIJRestoreArray_SeqAIJ" 2512 PetscErrorCode MatSeqAIJRestoreArray_SeqAIJ(Mat A,PetscScalar *array[]) 2513 { 2514 PetscFunctionBegin; 2515 PetscFunctionReturn(0); 2516 } 2517 2518 #undef __FUNCT__ 2519 #define __FUNCT__ "MatFDColoringApply_SeqAIJ" 2520 PetscErrorCode MatFDColoringApply_SeqAIJ(Mat J,MatFDColoring coloring,Vec x1,MatStructure *flag,void *sctx) 2521 { 2522 PetscErrorCode (*f)(void*,Vec,Vec,void*) = (PetscErrorCode (*)(void*,Vec,Vec,void *))coloring->f; 2523 PetscErrorCode ierr; 2524 PetscInt k,N,start,end,l,row,col,srow,**vscaleforrow; 2525 PetscScalar dx,*y,*xx,*w3_array; 2526 PetscScalar *vscale_array; 2527 PetscReal epsilon = coloring->error_rel,umin = coloring->umin; 2528 Vec w1,w2,w3; 2529 void *fctx = coloring->fctx; 2530 PetscBool flg = PETSC_FALSE; 2531 2532 PetscFunctionBegin; 2533 if (!coloring->w1) { 2534 ierr = VecDuplicate(x1,&coloring->w1);CHKERRQ(ierr); 2535 ierr = PetscLogObjectParent(coloring,coloring->w1);CHKERRQ(ierr); 2536 ierr = VecDuplicate(x1,&coloring->w2);CHKERRQ(ierr); 2537 ierr = PetscLogObjectParent(coloring,coloring->w2);CHKERRQ(ierr); 2538 ierr = VecDuplicate(x1,&coloring->w3);CHKERRQ(ierr); 2539 ierr = PetscLogObjectParent(coloring,coloring->w3);CHKERRQ(ierr); 2540 } 2541 w1 = coloring->w1; w2 = coloring->w2; w3 = coloring->w3; 2542 2543 ierr = MatSetUnfactored(J);CHKERRQ(ierr); 2544 ierr = PetscOptionsGetBool(((PetscObject)coloring)->prefix,"-mat_fd_coloring_dont_rezero",&flg,PETSC_NULL);CHKERRQ(ierr); 2545 if (flg) { 2546 ierr = PetscInfo(coloring,"Not calling MatZeroEntries()\n");CHKERRQ(ierr); 2547 } else { 2548 PetscBool assembled; 2549 ierr = MatAssembled(J,&assembled);CHKERRQ(ierr); 2550 if (assembled) { 2551 ierr = MatZeroEntries(J);CHKERRQ(ierr); 2552 } 2553 } 2554 2555 ierr = VecGetOwnershipRange(x1,&start,&end);CHKERRQ(ierr); 2556 ierr = VecGetSize(x1,&N);CHKERRQ(ierr); 2557 2558 if (!coloring->fset) { 2559 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2560 ierr = (*f)(sctx,x1,w1,fctx);CHKERRQ(ierr); 2561 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2562 } else { 2563 coloring->fset = PETSC_FALSE; 2564 } 2565 2566 /* 2567 Compute all the scale factors and share with other processors 2568 */ 2569 ierr = VecGetArray(x1,&xx);CHKERRQ(ierr);xx = xx - start; 2570 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr);vscale_array = vscale_array - start; 2571 for (k=0; k<coloring->ncolors; k++) { 2572 /* 2573 Loop over each column associated with color adding the 2574 perturbation to the vector w3. 2575 */ 2576 for (l=0; l<coloring->ncolumns[k]; l++) { 2577 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2578 dx = xx[col]; 2579 if (dx == 0.0) dx = 1.0; 2580 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2581 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2582 dx *= epsilon; 2583 vscale_array[col] = 1.0/dx; 2584 } 2585 } 2586 vscale_array = vscale_array + start;ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2587 ierr = VecGhostUpdateBegin(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2588 ierr = VecGhostUpdateEnd(coloring->vscale,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 2589 2590 /* ierr = VecView(coloring->vscale,PETSC_VIEWER_STDOUT_WORLD); 2591 ierr = VecView(x1,PETSC_VIEWER_STDOUT_WORLD);*/ 2592 2593 if (coloring->vscaleforrow) vscaleforrow = coloring->vscaleforrow; 2594 else vscaleforrow = coloring->columnsforrow; 2595 2596 ierr = VecGetArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2597 /* 2598 Loop over each color 2599 */ 2600 for (k=0; k<coloring->ncolors; k++) { 2601 coloring->currentcolor = k; 2602 ierr = VecCopy(x1,w3);CHKERRQ(ierr); 2603 ierr = VecGetArray(w3,&w3_array);CHKERRQ(ierr);w3_array = w3_array - start; 2604 /* 2605 Loop over each column associated with color adding the 2606 perturbation to the vector w3. 2607 */ 2608 for (l=0; l<coloring->ncolumns[k]; l++) { 2609 col = coloring->columns[k][l]; /* column of the matrix we are probing for */ 2610 dx = xx[col]; 2611 if (dx == 0.0) dx = 1.0; 2612 if (PetscAbsScalar(dx) < umin && PetscRealPart(dx) >= 0.0) dx = umin; 2613 else if (PetscRealPart(dx) < 0.0 && PetscAbsScalar(dx) < umin) dx = -umin; 2614 dx *= epsilon; 2615 if (!PetscAbsScalar(dx)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Computed 0 differencing parameter"); 2616 w3_array[col] += dx; 2617 } 2618 w3_array = w3_array + start; ierr = VecRestoreArray(w3,&w3_array);CHKERRQ(ierr); 2619 2620 /* 2621 Evaluate function at x1 + dx (here dx is a vector of perturbations) 2622 */ 2623 2624 ierr = PetscLogEventBegin(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2625 ierr = (*f)(sctx,w3,w2,fctx);CHKERRQ(ierr); 2626 ierr = PetscLogEventEnd(MAT_FDColoringFunction,0,0,0,0);CHKERRQ(ierr); 2627 ierr = VecAXPY(w2,-1.0,w1);CHKERRQ(ierr); 2628 2629 /* 2630 Loop over rows of vector, putting results into Jacobian matrix 2631 */ 2632 ierr = VecGetArray(w2,&y);CHKERRQ(ierr); 2633 for (l=0; l<coloring->nrows[k]; l++) { 2634 row = coloring->rows[k][l]; 2635 col = coloring->columnsforrow[k][l]; 2636 y[row] *= vscale_array[vscaleforrow[k][l]]; 2637 srow = row + start; 2638 ierr = MatSetValues_SeqAIJ(J,1,&srow,1,&col,y+row,INSERT_VALUES);CHKERRQ(ierr); 2639 } 2640 ierr = VecRestoreArray(w2,&y);CHKERRQ(ierr); 2641 } 2642 coloring->currentcolor = k; 2643 ierr = VecRestoreArray(coloring->vscale,&vscale_array);CHKERRQ(ierr); 2644 xx = xx + start; ierr = VecRestoreArray(x1,&xx);CHKERRQ(ierr); 2645 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2646 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 2647 PetscFunctionReturn(0); 2648 } 2649 2650 /* 2651 Computes the number of nonzeros per row needed for preallocation when X and Y 2652 have different nonzero structure. 2653 */ 2654 #undef __FUNCT__ 2655 #define __FUNCT__ "MatAXPYGetPreallocation_SeqAIJ" 2656 PetscErrorCode MatAXPYGetPreallocation_SeqAIJ(Mat Y,Mat X,PetscInt* nnz) 2657 { 2658 PetscInt i,m=Y->rmap->N; 2659 Mat_SeqAIJ *x = (Mat_SeqAIJ*)X->data; 2660 Mat_SeqAIJ *y = (Mat_SeqAIJ*)Y->data; 2661 const PetscInt *xi = x->i,*yi = y->i; 2662 2663 PetscFunctionBegin; 2664 /* Set the number of nonzeros in the new matrix */ 2665 for (i=0; i<m; i++) { 2666 PetscInt j,k,nzx = xi[i+1] - xi[i],nzy = yi[i+1] - yi[i]; 2667 const PetscInt *xj = x->j+xi[i],*yj = y->j+yi[i]; 2668 nnz[i] = 0; 2669 for (j=0,k=0; j<nzx; j++) { /* Point in X */ 2670 for (; k<nzy && yj[k]<xj[j]; k++) nnz[i]++; /* Catch up to X */ 2671 if (k<nzy && yj[k]==xj[j]) k++; /* Skip duplicate */ 2672 nnz[i]++; 2673 } 2674 for (; k<nzy; k++) nnz[i]++; 2675 } 2676 PetscFunctionReturn(0); 2677 } 2678 2679 #undef __FUNCT__ 2680 #define __FUNCT__ "MatAXPY_SeqAIJ" 2681 PetscErrorCode MatAXPY_SeqAIJ(Mat Y,PetscScalar a,Mat X,MatStructure str) 2682 { 2683 PetscErrorCode ierr; 2684 PetscInt i; 2685 Mat_SeqAIJ *x = (Mat_SeqAIJ *)X->data,*y = (Mat_SeqAIJ *)Y->data; 2686 PetscBLASInt one=1,bnz = PetscBLASIntCast(x->nz); 2687 2688 PetscFunctionBegin; 2689 if (str == SAME_NONZERO_PATTERN) { 2690 PetscScalar alpha = a; 2691 BLASaxpy_(&bnz,&alpha,x->a,&one,y->a,&one); 2692 y->idiagvalid = PETSC_FALSE; 2693 y->ibdiagvalid = PETSC_FALSE; 2694 } else if (str == SUBSET_NONZERO_PATTERN) { /* nonzeros of X is a subset of Y's */ 2695 if (y->xtoy && y->XtoY != X) { 2696 ierr = PetscFree(y->xtoy);CHKERRQ(ierr); 2697 ierr = MatDestroy(&y->XtoY);CHKERRQ(ierr); 2698 } 2699 if (!y->xtoy) { /* get xtoy */ 2700 ierr = MatAXPYGetxtoy_Private(X->rmap->n,x->i,x->j,PETSC_NULL, y->i,y->j,PETSC_NULL, &y->xtoy);CHKERRQ(ierr); 2701 y->XtoY = X; 2702 ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr); 2703 } 2704 for (i=0; i<x->nz; i++) y->a[y->xtoy[i]] += a*(x->a[i]); 2705 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); 2706 } else { 2707 Mat B; 2708 PetscInt *nnz; 2709 ierr = PetscMalloc(Y->rmap->N*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 2710 ierr = MatCreate(((PetscObject)Y)->comm,&B);CHKERRQ(ierr); 2711 ierr = PetscObjectSetName((PetscObject)B,((PetscObject)Y)->name);CHKERRQ(ierr); 2712 ierr = MatSetSizes(B,Y->rmap->n,Y->cmap->n,Y->rmap->N,Y->cmap->N);CHKERRQ(ierr); 2713 ierr = MatSetBlockSizes(B,Y->rmap->bs,Y->cmap->bs);CHKERRQ(ierr); 2714 ierr = MatSetType(B,(MatType) ((PetscObject)Y)->type_name);CHKERRQ(ierr); 2715 ierr = MatAXPYGetPreallocation_SeqAIJ(Y,X,nnz);CHKERRQ(ierr); 2716 ierr = MatSeqAIJSetPreallocation(B,0,nnz);CHKERRQ(ierr); 2717 ierr = MatAXPY_BasicWithPreallocation(B,Y,a,X,str);CHKERRQ(ierr); 2718 ierr = MatHeaderReplace(Y,B);CHKERRQ(ierr); 2719 ierr = PetscFree(nnz);CHKERRQ(ierr); 2720 } 2721 PetscFunctionReturn(0); 2722 } 2723 2724 #undef __FUNCT__ 2725 #define __FUNCT__ "MatConjugate_SeqAIJ" 2726 PetscErrorCode MatConjugate_SeqAIJ(Mat mat) 2727 { 2728 #if defined(PETSC_USE_COMPLEX) 2729 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 2730 PetscInt i,nz; 2731 PetscScalar *a; 2732 2733 PetscFunctionBegin; 2734 nz = aij->nz; 2735 a = aij->a; 2736 for (i=0; i<nz; i++) { 2737 a[i] = PetscConj(a[i]); 2738 } 2739 #else 2740 PetscFunctionBegin; 2741 #endif 2742 PetscFunctionReturn(0); 2743 } 2744 2745 #undef __FUNCT__ 2746 #define __FUNCT__ "MatGetRowMaxAbs_SeqAIJ" 2747 PetscErrorCode MatGetRowMaxAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2748 { 2749 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2750 PetscErrorCode ierr; 2751 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2752 PetscReal atmp; 2753 PetscScalar *x; 2754 MatScalar *aa; 2755 2756 PetscFunctionBegin; 2757 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2758 aa = a->a; 2759 ai = a->i; 2760 aj = a->j; 2761 2762 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2763 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2764 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2765 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2766 for (i=0; i<m; i++) { 2767 ncols = ai[1] - ai[0]; ai++; 2768 x[i] = 0.0; 2769 for (j=0; j<ncols; j++){ 2770 atmp = PetscAbsScalar(*aa); 2771 if (PetscAbsScalar(x[i]) < atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2772 aa++; aj++; 2773 } 2774 } 2775 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2776 PetscFunctionReturn(0); 2777 } 2778 2779 #undef __FUNCT__ 2780 #define __FUNCT__ "MatGetRowMax_SeqAIJ" 2781 PetscErrorCode MatGetRowMax_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2782 { 2783 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2784 PetscErrorCode ierr; 2785 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2786 PetscScalar *x; 2787 MatScalar *aa; 2788 2789 PetscFunctionBegin; 2790 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2791 aa = a->a; 2792 ai = a->i; 2793 aj = a->j; 2794 2795 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2796 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2797 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2798 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2799 for (i=0; i<m; i++) { 2800 ncols = ai[1] - ai[0]; ai++; 2801 if (ncols == A->cmap->n) { /* row is dense */ 2802 x[i] = *aa; if (idx) idx[i] = 0; 2803 } else { /* row is sparse so already KNOW maximum is 0.0 or higher */ 2804 x[i] = 0.0; 2805 if (idx) { 2806 idx[i] = 0; /* in case ncols is zero */ 2807 for (j=0;j<ncols;j++) { /* find first implicit 0.0 in the row */ 2808 if (aj[j] > j) { 2809 idx[i] = j; 2810 break; 2811 } 2812 } 2813 } 2814 } 2815 for (j=0; j<ncols; j++){ 2816 if (PetscRealPart(x[i]) < PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2817 aa++; aj++; 2818 } 2819 } 2820 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2821 PetscFunctionReturn(0); 2822 } 2823 2824 #undef __FUNCT__ 2825 #define __FUNCT__ "MatGetRowMinAbs_SeqAIJ" 2826 PetscErrorCode MatGetRowMinAbs_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2827 { 2828 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2829 PetscErrorCode ierr; 2830 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2831 PetscReal atmp; 2832 PetscScalar *x; 2833 MatScalar *aa; 2834 2835 PetscFunctionBegin; 2836 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2837 aa = a->a; 2838 ai = a->i; 2839 aj = a->j; 2840 2841 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2842 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2843 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2844 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); 2845 for (i=0; i<m; i++) { 2846 ncols = ai[1] - ai[0]; ai++; 2847 if (ncols) { 2848 /* Get first nonzero */ 2849 for (j = 0; j < ncols; j++) { 2850 atmp = PetscAbsScalar(aa[j]); 2851 if (atmp > 1.0e-12) {x[i] = atmp; if (idx) idx[i] = aj[j]; break;} 2852 } 2853 if (j == ncols) {x[i] = PetscAbsScalar(*aa); if (idx) idx[i] = *aj;} 2854 } else { 2855 x[i] = 0.0; if (idx) idx[i] = 0; 2856 } 2857 for (j = 0; j < ncols; j++) { 2858 atmp = PetscAbsScalar(*aa); 2859 if (atmp > 1.0e-12 && PetscAbsScalar(x[i]) > atmp) {x[i] = atmp; if (idx) idx[i] = *aj;} 2860 aa++; aj++; 2861 } 2862 } 2863 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2864 PetscFunctionReturn(0); 2865 } 2866 2867 #undef __FUNCT__ 2868 #define __FUNCT__ "MatGetRowMin_SeqAIJ" 2869 PetscErrorCode MatGetRowMin_SeqAIJ(Mat A,Vec v,PetscInt idx[]) 2870 { 2871 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 2872 PetscErrorCode ierr; 2873 PetscInt i,j,m = A->rmap->n,*ai,*aj,ncols,n; 2874 PetscScalar *x; 2875 MatScalar *aa; 2876 2877 PetscFunctionBegin; 2878 if (A->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 2879 aa = a->a; 2880 ai = a->i; 2881 aj = a->j; 2882 2883 ierr = VecSet(v,0.0);CHKERRQ(ierr); 2884 ierr = VecGetArray(v,&x);CHKERRQ(ierr); 2885 ierr = VecGetLocalSize(v,&n);CHKERRQ(ierr); 2886 if (n != A->rmap->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"Nonconforming matrix and vector"); 2887 for (i=0; i<m; i++) { 2888 ncols = ai[1] - ai[0]; ai++; 2889 if (ncols == A->cmap->n) { /* row is dense */ 2890 x[i] = *aa; if (idx) idx[i] = 0; 2891 } else { /* row is sparse so already KNOW minimum is 0.0 or lower */ 2892 x[i] = 0.0; 2893 if (idx) { /* find first implicit 0.0 in the row */ 2894 idx[i] = 0; /* in case ncols is zero */ 2895 for (j=0;j<ncols;j++) { 2896 if (aj[j] > j) { 2897 idx[i] = j; 2898 break; 2899 } 2900 } 2901 } 2902 } 2903 for (j=0; j<ncols; j++){ 2904 if (PetscRealPart(x[i]) > PetscRealPart(*aa)) {x[i] = *aa; if (idx) idx[i] = *aj;} 2905 aa++; aj++; 2906 } 2907 } 2908 ierr = VecRestoreArray(v,&x);CHKERRQ(ierr); 2909 PetscFunctionReturn(0); 2910 } 2911 2912 #include <petscblaslapack.h> 2913 #include <../src/mat/blockinvert.h> 2914 2915 #undef __FUNCT__ 2916 #define __FUNCT__ "MatInvertBlockDiagonal_SeqAIJ" 2917 PetscErrorCode MatInvertBlockDiagonal_SeqAIJ(Mat A,const PetscScalar **values) 2918 { 2919 Mat_SeqAIJ *a = (Mat_SeqAIJ*) A->data; 2920 PetscErrorCode ierr; 2921 PetscInt i,bs = A->rmap->bs,mbs = A->rmap->n/A->rmap->bs,ipvt[5],bs2 = bs*bs,*v_pivots,ij[7],*IJ,j; 2922 MatScalar *diag,work[25],*v_work; 2923 PetscReal shift = 0.0; 2924 2925 PetscFunctionBegin; 2926 if (a->ibdiagvalid) { 2927 if (values) *values = a->ibdiag; 2928 PetscFunctionReturn(0); 2929 } 2930 ierr = MatMarkDiagonal_SeqAIJ(A);CHKERRQ(ierr); 2931 if (!a->ibdiag) { 2932 ierr = PetscMalloc(bs2*mbs*sizeof(PetscScalar),&a->ibdiag);CHKERRQ(ierr); 2933 ierr = PetscLogObjectMemory(A,bs2*mbs*sizeof(PetscScalar));CHKERRQ(ierr); 2934 } 2935 diag = a->ibdiag; 2936 if (values) *values = a->ibdiag; 2937 /* factor and invert each block */ 2938 switch (bs){ 2939 case 1: 2940 for (i=0; i<mbs; i++) { 2941 ierr = MatGetValues(A,1,&i,1,&i,diag+i);CHKERRQ(ierr); 2942 diag[i] = (PetscScalar)1.0 / (diag[i] + shift); 2943 } 2944 break; 2945 case 2: 2946 for (i=0; i<mbs; i++) { 2947 ij[0] = 2*i; ij[1] = 2*i + 1; 2948 ierr = MatGetValues(A,2,ij,2,ij,diag);CHKERRQ(ierr); 2949 ierr = PetscKernel_A_gets_inverse_A_2(diag,shift);CHKERRQ(ierr); 2950 ierr = PetscKernel_A_gets_transpose_A_2(diag);CHKERRQ(ierr); 2951 diag += 4; 2952 } 2953 break; 2954 case 3: 2955 for (i=0; i<mbs; i++) { 2956 ij[0] = 3*i; ij[1] = 3*i + 1; ij[2] = 3*i + 2; 2957 ierr = MatGetValues(A,3,ij,3,ij,diag);CHKERRQ(ierr); 2958 ierr = PetscKernel_A_gets_inverse_A_3(diag,shift);CHKERRQ(ierr); 2959 ierr = PetscKernel_A_gets_transpose_A_3(diag);CHKERRQ(ierr); 2960 diag += 9; 2961 } 2962 break; 2963 case 4: 2964 for (i=0; i<mbs; i++) { 2965 ij[0] = 4*i; ij[1] = 4*i + 1; ij[2] = 4*i + 2; ij[3] = 4*i + 3; 2966 ierr = MatGetValues(A,4,ij,4,ij,diag);CHKERRQ(ierr); 2967 ierr = PetscKernel_A_gets_inverse_A_4(diag,shift);CHKERRQ(ierr); 2968 ierr = PetscKernel_A_gets_transpose_A_4(diag);CHKERRQ(ierr); 2969 diag += 16; 2970 } 2971 break; 2972 case 5: 2973 for (i=0; i<mbs; i++) { 2974 ij[0] = 5*i; ij[1] = 5*i + 1; ij[2] = 5*i + 2; ij[3] = 5*i + 3; ij[4] = 5*i + 4; 2975 ierr = MatGetValues(A,5,ij,5,ij,diag);CHKERRQ(ierr); 2976 ierr = PetscKernel_A_gets_inverse_A_5(diag,ipvt,work,shift);CHKERRQ(ierr); 2977 ierr = PetscKernel_A_gets_transpose_A_5(diag);CHKERRQ(ierr); 2978 diag += 25; 2979 } 2980 break; 2981 case 6: 2982 for (i=0; i<mbs; i++) { 2983 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; 2984 ierr = MatGetValues(A,6,ij,6,ij,diag);CHKERRQ(ierr); 2985 ierr = PetscKernel_A_gets_inverse_A_6(diag,shift);CHKERRQ(ierr); 2986 ierr = PetscKernel_A_gets_transpose_A_6(diag);CHKERRQ(ierr); 2987 diag += 36; 2988 } 2989 break; 2990 case 7: 2991 for (i=0; i<mbs; i++) { 2992 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; 2993 ierr = MatGetValues(A,7,ij,7,ij,diag);CHKERRQ(ierr); 2994 ierr = PetscKernel_A_gets_inverse_A_7(diag,shift);CHKERRQ(ierr); 2995 ierr = PetscKernel_A_gets_transpose_A_7(diag);CHKERRQ(ierr); 2996 diag += 49; 2997 } 2998 break; 2999 default: 3000 ierr = PetscMalloc3(bs,MatScalar,&v_work,bs,PetscInt,&v_pivots,bs,PetscInt,&IJ);CHKERRQ(ierr); 3001 for (i=0; i<mbs; i++) { 3002 for (j=0; j<bs; j++) { 3003 IJ[j] = bs*i + j; 3004 } 3005 ierr = MatGetValues(A,bs,IJ,bs,IJ,diag);CHKERRQ(ierr); 3006 ierr = PetscKernel_A_gets_inverse_A(bs,diag,v_pivots,v_work);CHKERRQ(ierr); 3007 ierr = PetscKernel_A_gets_transpose_A_N(diag,bs);CHKERRQ(ierr); 3008 diag += bs2; 3009 } 3010 ierr = PetscFree3(v_work,v_pivots,IJ);CHKERRQ(ierr); 3011 } 3012 a->ibdiagvalid = PETSC_TRUE; 3013 PetscFunctionReturn(0); 3014 } 3015 3016 #undef __FUNCT__ 3017 #define __FUNCT__ "MatSetRandom_SeqAIJ" 3018 static PetscErrorCode MatSetRandom_SeqAIJ(Mat x,PetscRandom rctx) 3019 { 3020 PetscErrorCode ierr; 3021 Mat_SeqAIJ *aij = (Mat_SeqAIJ*)x->data; 3022 PetscScalar a; 3023 PetscInt m,n,i,j,col; 3024 3025 PetscFunctionBegin; 3026 if (!x->assembled) { 3027 ierr = MatGetSize(x,&m,&n);CHKERRQ(ierr); 3028 for (i=0; i<m; i++) { 3029 for (j=0; j<aij->imax[i]; j++) { 3030 ierr = PetscRandomGetValue(rctx,&a);CHKERRQ(ierr); 3031 col = (PetscInt)(n*PetscRealPart(a)); 3032 ierr = MatSetValues(x,1,&i,1,&col,&a,ADD_VALUES);CHKERRQ(ierr); 3033 } 3034 } 3035 } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"Not yet coded"); 3036 ierr = MatAssemblyBegin(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3037 ierr = MatAssemblyEnd(x,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3038 PetscFunctionReturn(0); 3039 } 3040 3041 extern PetscErrorCode MatFDColoringApply_AIJ(Mat,MatFDColoring,Vec,MatStructure*,void*); 3042 /* -------------------------------------------------------------------*/ 3043 static struct _MatOps MatOps_Values = {MatSetValues_SeqAIJ, 3044 MatGetRow_SeqAIJ, 3045 MatRestoreRow_SeqAIJ, 3046 MatMult_SeqAIJ, 3047 /* 4*/ MatMultAdd_SeqAIJ, 3048 MatMultTranspose_SeqAIJ, 3049 MatMultTransposeAdd_SeqAIJ, 3050 0, 3051 0, 3052 0, 3053 /*10*/ 0, 3054 MatLUFactor_SeqAIJ, 3055 0, 3056 MatSOR_SeqAIJ, 3057 MatTranspose_SeqAIJ, 3058 /*15*/ MatGetInfo_SeqAIJ, 3059 MatEqual_SeqAIJ, 3060 MatGetDiagonal_SeqAIJ, 3061 MatDiagonalScale_SeqAIJ, 3062 MatNorm_SeqAIJ, 3063 /*20*/ 0, 3064 MatAssemblyEnd_SeqAIJ, 3065 MatSetOption_SeqAIJ, 3066 MatZeroEntries_SeqAIJ, 3067 /*24*/ MatZeroRows_SeqAIJ, 3068 0, 3069 0, 3070 0, 3071 0, 3072 /*29*/ MatSetUp_SeqAIJ, 3073 0, 3074 0, 3075 0, 3076 0, 3077 /*34*/ MatDuplicate_SeqAIJ, 3078 0, 3079 0, 3080 MatILUFactor_SeqAIJ, 3081 0, 3082 /*39*/ MatAXPY_SeqAIJ, 3083 MatGetSubMatrices_SeqAIJ, 3084 MatIncreaseOverlap_SeqAIJ, 3085 MatGetValues_SeqAIJ, 3086 MatCopy_SeqAIJ, 3087 /*44*/ MatGetRowMax_SeqAIJ, 3088 MatScale_SeqAIJ, 3089 0, 3090 MatDiagonalSet_SeqAIJ, 3091 MatZeroRowsColumns_SeqAIJ, 3092 /*49*/ MatSetRandom_SeqAIJ, 3093 MatGetRowIJ_SeqAIJ, 3094 MatRestoreRowIJ_SeqAIJ, 3095 MatGetColumnIJ_SeqAIJ, 3096 MatRestoreColumnIJ_SeqAIJ, 3097 /*54*/ MatFDColoringCreate_SeqAIJ, 3098 0, 3099 0, 3100 MatPermute_SeqAIJ, 3101 0, 3102 /*59*/ 0, 3103 MatDestroy_SeqAIJ, 3104 MatView_SeqAIJ, 3105 0, 3106 MatMatMatMult_SeqAIJ_SeqAIJ_SeqAIJ, 3107 /*64*/ MatMatMatMultSymbolic_SeqAIJ_SeqAIJ_SeqAIJ, 3108 MatMatMatMultNumeric_SeqAIJ_SeqAIJ_SeqAIJ, 3109 0, 3110 0, 3111 0, 3112 /*69*/ MatGetRowMaxAbs_SeqAIJ, 3113 MatGetRowMinAbs_SeqAIJ, 3114 0, 3115 MatSetColoring_SeqAIJ, 3116 0, 3117 /*74*/ MatSetValuesAdifor_SeqAIJ, 3118 MatFDColoringApply_AIJ, 3119 0, 3120 0, 3121 0, 3122 /*79*/ MatFindZeroDiagonals_SeqAIJ, 3123 0, 3124 0, 3125 0, 3126 MatLoad_SeqAIJ, 3127 /*84*/ MatIsSymmetric_SeqAIJ, 3128 MatIsHermitian_SeqAIJ, 3129 0, 3130 0, 3131 0, 3132 /*89*/ MatMatMult_SeqAIJ_SeqAIJ, 3133 MatMatMultSymbolic_SeqAIJ_SeqAIJ, 3134 MatMatMultNumeric_SeqAIJ_SeqAIJ, 3135 MatPtAP_SeqAIJ_SeqAIJ, 3136 MatPtAPSymbolic_SeqAIJ_SeqAIJ, 3137 /*94*/ MatPtAPNumeric_SeqAIJ_SeqAIJ, 3138 MatMatTransposeMult_SeqAIJ_SeqAIJ, 3139 MatMatTransposeMultSymbolic_SeqAIJ_SeqAIJ, 3140 MatMatTransposeMultNumeric_SeqAIJ_SeqAIJ, 3141 0, 3142 /*99*/ 0, 3143 0, 3144 0, 3145 MatConjugate_SeqAIJ, 3146 0, 3147 /*104*/MatSetValuesRow_SeqAIJ, 3148 MatRealPart_SeqAIJ, 3149 MatImaginaryPart_SeqAIJ, 3150 0, 3151 0, 3152 /*109*/MatMatSolve_SeqAIJ, 3153 0, 3154 MatGetRowMin_SeqAIJ, 3155 0, 3156 MatMissingDiagonal_SeqAIJ, 3157 /*114*/0, 3158 0, 3159 0, 3160 0, 3161 0, 3162 /*119*/0, 3163 0, 3164 0, 3165 0, 3166 MatGetMultiProcBlock_SeqAIJ, 3167 /*124*/MatFindNonzeroRows_SeqAIJ, 3168 MatGetColumnNorms_SeqAIJ, 3169 MatInvertBlockDiagonal_SeqAIJ, 3170 0, 3171 0, 3172 /*129*/0, 3173 MatTransposeMatMult_SeqAIJ_SeqAIJ, 3174 MatTransposeMatMultSymbolic_SeqAIJ_SeqAIJ, 3175 MatTransposeMatMultNumeric_SeqAIJ_SeqAIJ, 3176 MatTransposeColoringCreate_SeqAIJ, 3177 /*134*/MatTransColoringApplySpToDen_SeqAIJ, 3178 MatTransColoringApplyDenToSp_SeqAIJ, 3179 MatRARt_SeqAIJ_SeqAIJ, 3180 MatRARtSymbolic_SeqAIJ_SeqAIJ, 3181 MatRARtNumeric_SeqAIJ_SeqAIJ 3182 }; 3183 3184 EXTERN_C_BEGIN 3185 #undef __FUNCT__ 3186 #define __FUNCT__ "MatSeqAIJSetColumnIndices_SeqAIJ" 3187 PetscErrorCode MatSeqAIJSetColumnIndices_SeqAIJ(Mat mat,PetscInt *indices) 3188 { 3189 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3190 PetscInt i,nz,n; 3191 3192 PetscFunctionBegin; 3193 3194 nz = aij->maxnz; 3195 n = mat->rmap->n; 3196 for (i=0; i<nz; i++) { 3197 aij->j[i] = indices[i]; 3198 } 3199 aij->nz = nz; 3200 for (i=0; i<n; i++) { 3201 aij->ilen[i] = aij->imax[i]; 3202 } 3203 3204 PetscFunctionReturn(0); 3205 } 3206 EXTERN_C_END 3207 3208 #undef __FUNCT__ 3209 #define __FUNCT__ "MatSeqAIJSetColumnIndices" 3210 /*@ 3211 MatSeqAIJSetColumnIndices - Set the column indices for all the rows 3212 in the matrix. 3213 3214 Input Parameters: 3215 + mat - the SeqAIJ matrix 3216 - indices - the column indices 3217 3218 Level: advanced 3219 3220 Notes: 3221 This can be called if you have precomputed the nonzero structure of the 3222 matrix and want to provide it to the matrix object to improve the performance 3223 of the MatSetValues() operation. 3224 3225 You MUST have set the correct numbers of nonzeros per row in the call to 3226 MatCreateSeqAIJ(), and the columns indices MUST be sorted. 3227 3228 MUST be called before any calls to MatSetValues(); 3229 3230 The indices should start with zero, not one. 3231 3232 @*/ 3233 PetscErrorCode MatSeqAIJSetColumnIndices(Mat mat,PetscInt *indices) 3234 { 3235 PetscErrorCode ierr; 3236 3237 PetscFunctionBegin; 3238 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3239 PetscValidPointer(indices,2); 3240 ierr = PetscUseMethod(mat,"MatSeqAIJSetColumnIndices_C",(Mat,PetscInt *),(mat,indices));CHKERRQ(ierr); 3241 PetscFunctionReturn(0); 3242 } 3243 3244 /* ----------------------------------------------------------------------------------------*/ 3245 3246 EXTERN_C_BEGIN 3247 #undef __FUNCT__ 3248 #define __FUNCT__ "MatStoreValues_SeqAIJ" 3249 PetscErrorCode MatStoreValues_SeqAIJ(Mat mat) 3250 { 3251 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3252 PetscErrorCode ierr; 3253 size_t nz = aij->i[mat->rmap->n]; 3254 3255 PetscFunctionBegin; 3256 if (aij->nonew != 1) { 3257 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3258 } 3259 3260 /* allocate space for values if not already there */ 3261 if (!aij->saved_values) { 3262 ierr = PetscMalloc((nz+1)*sizeof(PetscScalar),&aij->saved_values);CHKERRQ(ierr); 3263 ierr = PetscLogObjectMemory(mat,(nz+1)*sizeof(PetscScalar));CHKERRQ(ierr); 3264 } 3265 3266 /* copy values over */ 3267 ierr = PetscMemcpy(aij->saved_values,aij->a,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3268 PetscFunctionReturn(0); 3269 } 3270 EXTERN_C_END 3271 3272 #undef __FUNCT__ 3273 #define __FUNCT__ "MatStoreValues" 3274 /*@ 3275 MatStoreValues - Stashes a copy of the matrix values; this allows, for 3276 example, reuse of the linear part of a Jacobian, while recomputing the 3277 nonlinear portion. 3278 3279 Collect on Mat 3280 3281 Input Parameters: 3282 . mat - the matrix (currently only AIJ matrices support this option) 3283 3284 Level: advanced 3285 3286 Common Usage, with SNESSolve(): 3287 $ Create Jacobian matrix 3288 $ Set linear terms into matrix 3289 $ Apply boundary conditions to matrix, at this time matrix must have 3290 $ final nonzero structure (i.e. setting the nonlinear terms and applying 3291 $ boundary conditions again will not change the nonzero structure 3292 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3293 $ ierr = MatStoreValues(mat); 3294 $ Call SNESSetJacobian() with matrix 3295 $ In your Jacobian routine 3296 $ ierr = MatRetrieveValues(mat); 3297 $ Set nonlinear terms in matrix 3298 3299 Common Usage without SNESSolve(), i.e. when you handle nonlinear solve yourself: 3300 $ // build linear portion of Jacobian 3301 $ ierr = MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); 3302 $ ierr = MatStoreValues(mat); 3303 $ loop over nonlinear iterations 3304 $ ierr = MatRetrieveValues(mat); 3305 $ // call MatSetValues(mat,...) to set nonliner portion of Jacobian 3306 $ // call MatAssemblyBegin/End() on matrix 3307 $ Solve linear system with Jacobian 3308 $ endloop 3309 3310 Notes: 3311 Matrix must already be assemblied before calling this routine 3312 Must set the matrix option MatSetOption(mat,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE); before 3313 calling this routine. 3314 3315 When this is called multiple times it overwrites the previous set of stored values 3316 and does not allocated additional space. 3317 3318 .seealso: MatRetrieveValues() 3319 3320 @*/ 3321 PetscErrorCode MatStoreValues(Mat mat) 3322 { 3323 PetscErrorCode ierr; 3324 3325 PetscFunctionBegin; 3326 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3327 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3328 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3329 ierr = PetscUseMethod(mat,"MatStoreValues_C",(Mat),(mat));CHKERRQ(ierr); 3330 PetscFunctionReturn(0); 3331 } 3332 3333 EXTERN_C_BEGIN 3334 #undef __FUNCT__ 3335 #define __FUNCT__ "MatRetrieveValues_SeqAIJ" 3336 PetscErrorCode MatRetrieveValues_SeqAIJ(Mat mat) 3337 { 3338 Mat_SeqAIJ *aij = (Mat_SeqAIJ *)mat->data; 3339 PetscErrorCode ierr; 3340 PetscInt nz = aij->i[mat->rmap->n]; 3341 3342 PetscFunctionBegin; 3343 if (aij->nonew != 1) { 3344 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatSetOption(A,MAT_NEW_NONZERO_LOCATIONS,PETSC_FALSE);first"); 3345 } 3346 if (!aij->saved_values) { 3347 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ORDER,"Must call MatStoreValues(A);first"); 3348 } 3349 /* copy values over */ 3350 ierr = PetscMemcpy(aij->a,aij->saved_values,nz*sizeof(PetscScalar));CHKERRQ(ierr); 3351 PetscFunctionReturn(0); 3352 } 3353 EXTERN_C_END 3354 3355 #undef __FUNCT__ 3356 #define __FUNCT__ "MatRetrieveValues" 3357 /*@ 3358 MatRetrieveValues - Retrieves the copy of the matrix values; this allows, for 3359 example, reuse of the linear part of a Jacobian, while recomputing the 3360 nonlinear portion. 3361 3362 Collect on Mat 3363 3364 Input Parameters: 3365 . mat - the matrix (currently on AIJ matrices support this option) 3366 3367 Level: advanced 3368 3369 .seealso: MatStoreValues() 3370 3371 @*/ 3372 PetscErrorCode MatRetrieveValues(Mat mat) 3373 { 3374 PetscErrorCode ierr; 3375 3376 PetscFunctionBegin; 3377 PetscValidHeaderSpecific(mat,MAT_CLASSID,1); 3378 if (!mat->assembled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for unassembled matrix"); 3379 if (mat->factortype) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Not for factored matrix"); 3380 ierr = PetscUseMethod(mat,"MatRetrieveValues_C",(Mat),(mat));CHKERRQ(ierr); 3381 PetscFunctionReturn(0); 3382 } 3383 3384 3385 /* --------------------------------------------------------------------------------*/ 3386 #undef __FUNCT__ 3387 #define __FUNCT__ "MatCreateSeqAIJ" 3388 /*@C 3389 MatCreateSeqAIJ - Creates a sparse matrix in AIJ (compressed row) format 3390 (the default parallel PETSc format). For good matrix assembly performance 3391 the user should preallocate the matrix storage by setting the parameter nz 3392 (or the array nnz). By setting these parameters accurately, performance 3393 during matrix assembly can be increased by more than a factor of 50. 3394 3395 Collective on MPI_Comm 3396 3397 Input Parameters: 3398 + comm - MPI communicator, set to PETSC_COMM_SELF 3399 . m - number of rows 3400 . n - number of columns 3401 . nz - number of nonzeros per row (same for all rows) 3402 - nnz - array containing the number of nonzeros in the various rows 3403 (possibly different for each row) or PETSC_NULL 3404 3405 Output Parameter: 3406 . A - the matrix 3407 3408 It is recommended that one use the MatCreate(), MatSetType() and/or MatSetFromOptions(), 3409 MatXXXXSetPreallocation() paradgm instead of this routine directly. 3410 [MatXXXXSetPreallocation() is, for example, MatSeqAIJSetPreallocation] 3411 3412 Notes: 3413 If nnz is given then nz is ignored 3414 3415 The AIJ format (also called the Yale sparse matrix format or 3416 compressed row storage), is fully compatible with standard Fortran 77 3417 storage. That is, the stored row and column indices can begin at 3418 either one (as in Fortran) or zero. See the users' manual for details. 3419 3420 Specify the preallocated storage with either nz or nnz (not both). 3421 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3422 allocation. For large problems you MUST preallocate memory or you 3423 will get TERRIBLE performance, see the users' manual chapter on matrices. 3424 3425 By default, this format uses inodes (identical nodes) when possible, to 3426 improve numerical efficiency of matrix-vector products and solves. We 3427 search for consecutive rows with the same nonzero structure, thereby 3428 reusing matrix information to achieve increased efficiency. 3429 3430 Options Database Keys: 3431 + -mat_no_inode - Do not use inodes 3432 - -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3433 3434 Level: intermediate 3435 3436 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays() 3437 3438 @*/ 3439 PetscErrorCode MatCreateSeqAIJ(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt nz,const PetscInt nnz[],Mat *A) 3440 { 3441 PetscErrorCode ierr; 3442 3443 PetscFunctionBegin; 3444 ierr = MatCreate(comm,A);CHKERRQ(ierr); 3445 ierr = MatSetSizes(*A,m,n,m,n);CHKERRQ(ierr); 3446 ierr = MatSetType(*A,MATSEQAIJ);CHKERRQ(ierr); 3447 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*A,nz,nnz);CHKERRQ(ierr); 3448 PetscFunctionReturn(0); 3449 } 3450 3451 #undef __FUNCT__ 3452 #define __FUNCT__ "MatSeqAIJSetPreallocation" 3453 /*@C 3454 MatSeqAIJSetPreallocation - For good matrix assembly performance 3455 the user should preallocate the matrix storage by setting the parameter nz 3456 (or the array nnz). By setting these parameters accurately, performance 3457 during matrix assembly can be increased by more than a factor of 50. 3458 3459 Collective on MPI_Comm 3460 3461 Input Parameters: 3462 + B - The matrix-free 3463 . nz - number of nonzeros per row (same for all rows) 3464 - nnz - array containing the number of nonzeros in the various rows 3465 (possibly different for each row) or PETSC_NULL 3466 3467 Notes: 3468 If nnz is given then nz is ignored 3469 3470 The AIJ format (also called the Yale sparse matrix format or 3471 compressed row storage), is fully compatible with standard Fortran 77 3472 storage. That is, the stored row and column indices can begin at 3473 either one (as in Fortran) or zero. See the users' manual for details. 3474 3475 Specify the preallocated storage with either nz or nnz (not both). 3476 Set nz=PETSC_DEFAULT and nnz=PETSC_NULL for PETSc to control dynamic memory 3477 allocation. For large problems you MUST preallocate memory or you 3478 will get TERRIBLE performance, see the users' manual chapter on matrices. 3479 3480 You can call MatGetInfo() to get information on how effective the preallocation was; 3481 for example the fields mallocs,nz_allocated,nz_used,nz_unneeded; 3482 You can also run with the option -info and look for messages with the string 3483 malloc in them to see if additional memory allocation was needed. 3484 3485 Developers: Use nz of MAT_SKIP_ALLOCATION to not allocate any space for the matrix 3486 entries or columns indices 3487 3488 By default, this format uses inodes (identical nodes) when possible, to 3489 improve numerical efficiency of matrix-vector products and solves. We 3490 search for consecutive rows with the same nonzero structure, thereby 3491 reusing matrix information to achieve increased efficiency. 3492 3493 Options Database Keys: 3494 + -mat_no_inode - Do not use inodes 3495 . -mat_inode_limit <limit> - Sets inode limit (max limit=5) 3496 - -mat_aij_oneindex - Internally use indexing starting at 1 3497 rather than 0. Note that when calling MatSetValues(), 3498 the user still MUST index entries starting at 0! 3499 3500 Level: intermediate 3501 3502 .seealso: MatCreate(), MatCreateAIJ(), MatSetValues(), MatSeqAIJSetColumnIndices(), MatCreateSeqAIJWithArrays(), MatGetInfo() 3503 3504 @*/ 3505 PetscErrorCode MatSeqAIJSetPreallocation(Mat B,PetscInt nz,const PetscInt nnz[]) 3506 { 3507 PetscErrorCode ierr; 3508 3509 PetscFunctionBegin; 3510 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3511 PetscValidType(B,1); 3512 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocation_C",(Mat,PetscInt,const PetscInt[]),(B,nz,nnz));CHKERRQ(ierr); 3513 PetscFunctionReturn(0); 3514 } 3515 3516 EXTERN_C_BEGIN 3517 #undef __FUNCT__ 3518 #define __FUNCT__ "MatSeqAIJSetPreallocation_SeqAIJ" 3519 PetscErrorCode MatSeqAIJSetPreallocation_SeqAIJ(Mat B,PetscInt nz,const PetscInt *nnz) 3520 { 3521 Mat_SeqAIJ *b; 3522 PetscBool skipallocation = PETSC_FALSE,realalloc = PETSC_FALSE; 3523 PetscErrorCode ierr; 3524 PetscInt i; 3525 3526 PetscFunctionBegin; 3527 if (nz >= 0 || nnz) realalloc = PETSC_TRUE; 3528 if (nz == MAT_SKIP_ALLOCATION) { 3529 skipallocation = PETSC_TRUE; 3530 nz = 0; 3531 } 3532 3533 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3534 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3535 3536 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 5; 3537 if (nz < 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"nz cannot be less than 0: value %d",nz); 3538 if (nnz) { 3539 for (i=0; i<B->rmap->n; i++) { 3540 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]); 3541 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); 3542 } 3543 } 3544 3545 B->preallocated = PETSC_TRUE; 3546 b = (Mat_SeqAIJ*)B->data; 3547 3548 if (!skipallocation) { 3549 if (!b->imax) { 3550 ierr = PetscMalloc2(B->rmap->n,PetscInt,&b->imax,B->rmap->n,PetscInt,&b->ilen);CHKERRQ(ierr); 3551 ierr = PetscLogObjectMemory(B,2*B->rmap->n*sizeof(PetscInt));CHKERRQ(ierr); 3552 } 3553 if (!nnz) { 3554 if (nz == PETSC_DEFAULT || nz == PETSC_DECIDE) nz = 10; 3555 else if (nz < 0) nz = 1; 3556 for (i=0; i<B->rmap->n; i++) b->imax[i] = nz; 3557 nz = nz*B->rmap->n; 3558 } else { 3559 nz = 0; 3560 for (i=0; i<B->rmap->n; i++) {b->imax[i] = nnz[i]; nz += nnz[i];} 3561 } 3562 /* b->ilen will count nonzeros in each row so far. */ 3563 for (i=0; i<B->rmap->n; i++) { b->ilen[i] = 0; } 3564 3565 /* allocate the matrix space */ 3566 ierr = MatSeqXAIJFreeAIJ(B,&b->a,&b->j,&b->i);CHKERRQ(ierr); 3567 ierr = PetscMalloc3(nz,PetscScalar,&b->a,nz,PetscInt,&b->j,B->rmap->n+1,PetscInt,&b->i);CHKERRQ(ierr); 3568 ierr = PetscLogObjectMemory(B,(B->rmap->n+1)*sizeof(PetscInt)+nz*(sizeof(PetscScalar)+sizeof(PetscInt)));CHKERRQ(ierr); 3569 b->i[0] = 0; 3570 for (i=1; i<B->rmap->n+1; i++) { 3571 b->i[i] = b->i[i-1] + b->imax[i-1]; 3572 } 3573 b->singlemalloc = PETSC_TRUE; 3574 b->free_a = PETSC_TRUE; 3575 b->free_ij = PETSC_TRUE; 3576 #if defined(PETSC_THREADCOMM_ACTIVE) 3577 ierr = MatZeroEntries_SeqAIJ(B);CHKERRQ(ierr); 3578 #endif 3579 } else { 3580 b->free_a = PETSC_FALSE; 3581 b->free_ij = PETSC_FALSE; 3582 } 3583 3584 b->nz = 0; 3585 b->maxnz = nz; 3586 B->info.nz_unneeded = (double)b->maxnz; 3587 if (realalloc) {ierr = MatSetOption(B,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr);} 3588 PetscFunctionReturn(0); 3589 } 3590 EXTERN_C_END 3591 3592 #undef __FUNCT__ 3593 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR" 3594 /*@ 3595 MatSeqAIJSetPreallocationCSR - Allocates memory for a sparse sequential matrix in AIJ format. 3596 3597 Input Parameters: 3598 + B - the matrix 3599 . i - the indices into j for the start of each row (starts with zero) 3600 . j - the column indices for each row (starts with zero) these must be sorted for each row 3601 - v - optional values in the matrix 3602 3603 Level: developer 3604 3605 The i,j,v values are COPIED with this routine; to avoid the copy use MatCreateSeqAIJWithArrays() 3606 3607 .keywords: matrix, aij, compressed row, sparse, sequential 3608 3609 .seealso: MatCreate(), MatCreateSeqAIJ(), MatSetValues(), MatSeqAIJSetPreallocation(), MatCreateSeqAIJ(), SeqAIJ 3610 @*/ 3611 PetscErrorCode MatSeqAIJSetPreallocationCSR(Mat B,const PetscInt i[],const PetscInt j[],const PetscScalar v[]) 3612 { 3613 PetscErrorCode ierr; 3614 3615 PetscFunctionBegin; 3616 PetscValidHeaderSpecific(B,MAT_CLASSID,1); 3617 PetscValidType(B,1); 3618 ierr = PetscTryMethod(B,"MatSeqAIJSetPreallocationCSR_C",(Mat,const PetscInt[],const PetscInt[],const PetscScalar[]),(B,i,j,v));CHKERRQ(ierr); 3619 PetscFunctionReturn(0); 3620 } 3621 3622 EXTERN_C_BEGIN 3623 #undef __FUNCT__ 3624 #define __FUNCT__ "MatSeqAIJSetPreallocationCSR_SeqAIJ" 3625 PetscErrorCode MatSeqAIJSetPreallocationCSR_SeqAIJ(Mat B,const PetscInt Ii[],const PetscInt J[],const PetscScalar v[]) 3626 { 3627 PetscInt i; 3628 PetscInt m,n; 3629 PetscInt nz; 3630 PetscInt *nnz, nz_max = 0; 3631 PetscScalar *values; 3632 PetscErrorCode ierr; 3633 3634 PetscFunctionBegin; 3635 if (Ii[0]) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Ii[0] must be 0 it is %D", Ii[0]); 3636 3637 ierr = PetscLayoutSetUp(B->rmap);CHKERRQ(ierr); 3638 ierr = PetscLayoutSetUp(B->cmap);CHKERRQ(ierr); 3639 3640 ierr = MatGetSize(B, &m, &n);CHKERRQ(ierr); 3641 ierr = PetscMalloc((m+1) * sizeof(PetscInt), &nnz);CHKERRQ(ierr); 3642 for (i = 0; i < m; i++) { 3643 nz = Ii[i+1]- Ii[i]; 3644 nz_max = PetscMax(nz_max, nz); 3645 if (nz < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE, "Local row %D has a negative number of columns %D", i, nnz); 3646 nnz[i] = nz; 3647 } 3648 ierr = MatSeqAIJSetPreallocation(B, 0, nnz);CHKERRQ(ierr); 3649 ierr = PetscFree(nnz);CHKERRQ(ierr); 3650 3651 if (v) { 3652 values = (PetscScalar*) v; 3653 } else { 3654 ierr = PetscMalloc(nz_max*sizeof(PetscScalar), &values);CHKERRQ(ierr); 3655 ierr = PetscMemzero(values, nz_max*sizeof(PetscScalar));CHKERRQ(ierr); 3656 } 3657 3658 for (i = 0; i < m; i++) { 3659 nz = Ii[i+1] - Ii[i]; 3660 ierr = MatSetValues_SeqAIJ(B, 1, &i, nz, J+Ii[i], values + (v ? Ii[i] : 0), INSERT_VALUES);CHKERRQ(ierr); 3661 } 3662 3663 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3664 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 3665 3666 if (!v) { 3667 ierr = PetscFree(values);CHKERRQ(ierr); 3668 } 3669 ierr = MatSetOption(B,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 3670 PetscFunctionReturn(0); 3671 } 3672 EXTERN_C_END 3673 3674 #include <../src/mat/impls/dense/seq/dense.h> 3675 #include <petsc-private/petscaxpy.h> 3676 3677 #undef __FUNCT__ 3678 #define __FUNCT__ "MatMatMultNumeric_SeqDense_SeqAIJ" 3679 /* 3680 Computes (B'*A')' since computing B*A directly is untenable 3681 3682 n p p 3683 ( ) ( ) ( ) 3684 m ( A ) * n ( B ) = m ( C ) 3685 ( ) ( ) ( ) 3686 3687 */ 3688 PetscErrorCode MatMatMultNumeric_SeqDense_SeqAIJ(Mat A,Mat B,Mat C) 3689 { 3690 PetscErrorCode ierr; 3691 Mat_SeqDense *sub_a = (Mat_SeqDense*)A->data; 3692 Mat_SeqAIJ *sub_b = (Mat_SeqAIJ*)B->data; 3693 Mat_SeqDense *sub_c = (Mat_SeqDense*)C->data; 3694 PetscInt i,n,m,q,p; 3695 const PetscInt *ii,*idx; 3696 const PetscScalar *b,*a,*a_q; 3697 PetscScalar *c,*c_q; 3698 3699 PetscFunctionBegin; 3700 m = A->rmap->n; 3701 n = A->cmap->n; 3702 p = B->cmap->n; 3703 a = sub_a->v; 3704 b = sub_b->a; 3705 c = sub_c->v; 3706 ierr = PetscMemzero(c,m*p*sizeof(PetscScalar));CHKERRQ(ierr); 3707 3708 ii = sub_b->i; 3709 idx = sub_b->j; 3710 for (i=0; i<n; i++) { 3711 q = ii[i+1] - ii[i]; 3712 while (q-->0) { 3713 c_q = c + m*(*idx); 3714 a_q = a + m*i; 3715 PetscAXPY(c_q,*b,a_q,m); 3716 idx++; 3717 b++; 3718 } 3719 } 3720 PetscFunctionReturn(0); 3721 } 3722 3723 #undef __FUNCT__ 3724 #define __FUNCT__ "MatMatMultSymbolic_SeqDense_SeqAIJ" 3725 PetscErrorCode MatMatMultSymbolic_SeqDense_SeqAIJ(Mat A,Mat B,PetscReal fill,Mat *C) 3726 { 3727 PetscErrorCode ierr; 3728 PetscInt m=A->rmap->n,n=B->cmap->n; 3729 Mat Cmat; 3730 3731 PetscFunctionBegin; 3732 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); 3733 ierr = MatCreate(((PetscObject)A)->comm,&Cmat);CHKERRQ(ierr); 3734 ierr = MatSetSizes(Cmat,m,n,m,n);CHKERRQ(ierr); 3735 ierr = MatSetBlockSizes(Cmat,A->rmap->bs,B->cmap->bs);CHKERRQ(ierr); 3736 ierr = MatSetType(Cmat,MATSEQDENSE);CHKERRQ(ierr); 3737 ierr = MatSeqDenseSetPreallocation(Cmat,PETSC_NULL);CHKERRQ(ierr); 3738 3739 Cmat->ops->matmultnumeric = MatMatMultNumeric_SeqDense_SeqAIJ; 3740 *C = Cmat; 3741 PetscFunctionReturn(0); 3742 } 3743 3744 /* ----------------------------------------------------------------*/ 3745 #undef __FUNCT__ 3746 #define __FUNCT__ "MatMatMult_SeqDense_SeqAIJ" 3747 PetscErrorCode MatMatMult_SeqDense_SeqAIJ(Mat A,Mat B,MatReuse scall,PetscReal fill,Mat *C) 3748 { 3749 PetscErrorCode ierr; 3750 3751 PetscFunctionBegin; 3752 if (scall == MAT_INITIAL_MATRIX){ 3753 ierr = MatMatMultSymbolic_SeqDense_SeqAIJ(A,B,fill,C);CHKERRQ(ierr); 3754 } 3755 ierr = MatMatMultNumeric_SeqDense_SeqAIJ(A,B,*C);CHKERRQ(ierr); 3756 PetscFunctionReturn(0); 3757 } 3758 3759 3760 /*MC 3761 MATSEQAIJ - MATSEQAIJ = "seqaij" - A matrix type to be used for sequential sparse matrices, 3762 based on compressed sparse row format. 3763 3764 Options Database Keys: 3765 . -mat_type seqaij - sets the matrix type to "seqaij" during a call to MatSetFromOptions() 3766 3767 Level: beginner 3768 3769 .seealso: MatCreateSeqAIJ(), MatSetFromOptions(), MatSetType(), MatCreate(), MatType 3770 M*/ 3771 3772 /*MC 3773 MATAIJ - MATAIJ = "aij" - A matrix type to be used for sparse matrices. 3774 3775 This matrix type is identical to MATSEQAIJ when constructed with a single process communicator, 3776 and MATMPIAIJ otherwise. As a result, for single process communicators, 3777 MatSeqAIJSetPreallocation is supported, and similarly MatMPIAIJSetPreallocation is supported 3778 for communicators controlling multiple processes. It is recommended that you call both of 3779 the above preallocation routines for simplicity. 3780 3781 Options Database Keys: 3782 . -mat_type aij - sets the matrix type to "aij" during a call to MatSetFromOptions() 3783 3784 Developer Notes: Subclasses include MATAIJCUSP, MATAIJPERM, MATAIJCRL, and also automatically switches over to use inodes when 3785 enough exist. 3786 3787 Level: beginner 3788 3789 .seealso: MatCreateAIJ(), MatCreateSeqAIJ(), MATSEQAIJ,MATMPIAIJ 3790 M*/ 3791 3792 /*MC 3793 MATAIJCRL - MATAIJCRL = "aijcrl" - A matrix type to be used for sparse matrices. 3794 3795 This matrix type is identical to MATSEQAIJCRL when constructed with a single process communicator, 3796 and MATMPIAIJCRL otherwise. As a result, for single process communicators, 3797 MatSeqAIJSetPreallocation() is supported, and similarly MatMPIAIJSetPreallocation() is supported 3798 for communicators controlling multiple processes. It is recommended that you call both of 3799 the above preallocation routines for simplicity. 3800 3801 Options Database Keys: 3802 . -mat_type aijcrl - sets the matrix type to "aijcrl" during a call to MatSetFromOptions() 3803 3804 Level: beginner 3805 3806 .seealso: MatCreateMPIAIJCRL,MATSEQAIJCRL,MATMPIAIJCRL, MATSEQAIJCRL, MATMPIAIJCRL 3807 M*/ 3808 3809 EXTERN_C_BEGIN 3810 #if defined(PETSC_HAVE_PASTIX) 3811 extern PetscErrorCode MatGetFactor_seqaij_pastix(Mat,MatFactorType,Mat*); 3812 #endif 3813 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3814 extern PetscErrorCode MatGetFactor_seqaij_essl(Mat,MatFactorType,Mat *); 3815 #endif 3816 extern PetscErrorCode MatConvert_SeqAIJ_SeqAIJCRL(Mat,MatType,MatReuse,Mat*); 3817 extern PetscErrorCode MatGetFactor_seqaij_petsc(Mat,MatFactorType,Mat*); 3818 extern PetscErrorCode MatGetFactor_seqaij_bas(Mat,MatFactorType,Mat*); 3819 extern PetscErrorCode MatGetFactorAvailable_seqaij_petsc(Mat,MatFactorType,PetscBool *); 3820 #if defined(PETSC_HAVE_MUMPS) 3821 extern PetscErrorCode MatGetFactor_aij_mumps(Mat,MatFactorType,Mat*); 3822 #endif 3823 #if defined(PETSC_HAVE_SUPERLU) 3824 extern PetscErrorCode MatGetFactor_seqaij_superlu(Mat,MatFactorType,Mat*); 3825 #endif 3826 #if defined(PETSC_HAVE_SUPERLU_DIST) 3827 extern PetscErrorCode MatGetFactor_seqaij_superlu_dist(Mat,MatFactorType,Mat*); 3828 #endif 3829 #if defined(PETSC_HAVE_UMFPACK) 3830 extern PetscErrorCode MatGetFactor_seqaij_umfpack(Mat,MatFactorType,Mat*); 3831 #endif 3832 #if defined(PETSC_HAVE_CHOLMOD) 3833 extern PetscErrorCode MatGetFactor_seqaij_cholmod(Mat,MatFactorType,Mat*); 3834 #endif 3835 #if defined(PETSC_HAVE_LUSOL) 3836 extern PetscErrorCode MatGetFactor_seqaij_lusol(Mat,MatFactorType,Mat*); 3837 #endif 3838 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3839 extern PetscErrorCode MatGetFactor_seqaij_matlab(Mat,MatFactorType,Mat*); 3840 extern PetscErrorCode MatlabEnginePut_SeqAIJ(PetscObject,void*); 3841 extern PetscErrorCode MatlabEngineGet_SeqAIJ(PetscObject,void*); 3842 #endif 3843 #if defined(PETSC_HAVE_CLIQUE) 3844 extern PetscErrorCode MatGetFactor_aij_clique(Mat,MatFactorType,Mat*); 3845 #endif 3846 EXTERN_C_END 3847 3848 3849 #undef __FUNCT__ 3850 #define __FUNCT__ "MatSeqAIJGetArray" 3851 /*@C 3852 MatSeqAIJGetArray - gives access to the array where the data for a SeqSeqAIJ matrix is stored 3853 3854 Not Collective 3855 3856 Input Parameter: 3857 . mat - a MATSEQDENSE matrix 3858 3859 Output Parameter: 3860 . array - pointer to the data 3861 3862 Level: intermediate 3863 3864 .seealso: MatSeqAIJRestoreArray() 3865 @*/ 3866 PetscErrorCode MatSeqAIJGetArray(Mat A,PetscScalar **array) 3867 { 3868 PetscErrorCode ierr; 3869 3870 PetscFunctionBegin; 3871 ierr = PetscUseMethod(A,"MatSeqAIJGetArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3872 PetscFunctionReturn(0); 3873 } 3874 3875 #undef __FUNCT__ 3876 #define __FUNCT__ "MatSeqAIJRestoreArray" 3877 /*@C 3878 MatSeqAIJRestoreArray - returns access to the array where the data for a SeqSeqAIJ matrix is stored obtained by MatSeqAIJGetArray() 3879 3880 Not Collective 3881 3882 Input Parameters: 3883 . mat - a MATSEQDENSE matrix 3884 . array - pointer to the data 3885 3886 Level: intermediate 3887 3888 .seealso: MatSeqAIJGetArray() 3889 @*/ 3890 PetscErrorCode MatSeqAIJRestoreArray(Mat A,PetscScalar **array) 3891 { 3892 PetscErrorCode ierr; 3893 3894 PetscFunctionBegin; 3895 ierr = PetscUseMethod(A,"MatSeqAIJRestoreArray_C",(Mat,PetscScalar**),(A,array));CHKERRQ(ierr); 3896 PetscFunctionReturn(0); 3897 } 3898 3899 EXTERN_C_BEGIN 3900 #undef __FUNCT__ 3901 #define __FUNCT__ "MatCreate_SeqAIJ" 3902 PetscErrorCode MatCreate_SeqAIJ(Mat B) 3903 { 3904 Mat_SeqAIJ *b; 3905 PetscErrorCode ierr; 3906 PetscMPIInt size; 3907 3908 PetscFunctionBegin; 3909 ierr = MPI_Comm_size(((PetscObject)B)->comm,&size);CHKERRQ(ierr); 3910 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Comm must be of size 1"); 3911 3912 ierr = PetscNewLog(B,Mat_SeqAIJ,&b);CHKERRQ(ierr); 3913 B->data = (void*)b; 3914 ierr = PetscMemcpy(B->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr); 3915 b->row = 0; 3916 b->col = 0; 3917 b->icol = 0; 3918 b->reallocs = 0; 3919 b->ignorezeroentries = PETSC_FALSE; 3920 b->roworiented = PETSC_TRUE; 3921 b->nonew = 0; 3922 b->diag = 0; 3923 b->solve_work = 0; 3924 B->spptr = 0; 3925 b->saved_values = 0; 3926 b->idiag = 0; 3927 b->mdiag = 0; 3928 b->ssor_work = 0; 3929 b->omega = 1.0; 3930 b->fshift = 0.0; 3931 b->idiagvalid = PETSC_FALSE; 3932 b->ibdiagvalid = PETSC_FALSE; 3933 b->keepnonzeropattern = PETSC_FALSE; 3934 b->xtoy = 0; 3935 b->XtoY = 0; 3936 B->same_nonzero = PETSC_FALSE; 3937 3938 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3939 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJGetArray_C","MatSeqAIJGetArray_SeqAIJ",MatSeqAIJGetArray_SeqAIJ);CHKERRQ(ierr); 3940 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJRestoreArray_C","MatSeqAIJRestoreArray_SeqAIJ",MatSeqAIJRestoreArray_SeqAIJ);CHKERRQ(ierr); 3941 3942 #if defined(PETSC_HAVE_MATLAB_ENGINE) 3943 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_matlab_C","MatGetFactor_seqaij_matlab",MatGetFactor_seqaij_matlab);CHKERRQ(ierr); 3944 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEnginePut_C","MatlabEnginePut_SeqAIJ",MatlabEnginePut_SeqAIJ);CHKERRQ(ierr); 3945 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"PetscMatlabEngineGet_C","MatlabEngineGet_SeqAIJ",MatlabEngineGet_SeqAIJ);CHKERRQ(ierr); 3946 #endif 3947 #if defined(PETSC_HAVE_PASTIX) 3948 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_pastix_C","MatGetFactor_seqaij_pastix",MatGetFactor_seqaij_pastix);CHKERRQ(ierr); 3949 #endif 3950 #if defined(PETSC_HAVE_ESSL) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_REAL_SINGLE) && !defined(PETSC_USE_REAL___FLOAT128) 3951 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_essl_C","MatGetFactor_seqaij_essl",MatGetFactor_seqaij_essl);CHKERRQ(ierr); 3952 #endif 3953 #if defined(PETSC_HAVE_SUPERLU) 3954 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_C","MatGetFactor_seqaij_superlu",MatGetFactor_seqaij_superlu);CHKERRQ(ierr); 3955 #endif 3956 #if defined(PETSC_HAVE_SUPERLU_DIST) 3957 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_superlu_dist_C","MatGetFactor_seqaij_superlu_dist",MatGetFactor_seqaij_superlu_dist);CHKERRQ(ierr); 3958 #endif 3959 #if defined(PETSC_HAVE_MUMPS) 3960 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_mumps_C","MatGetFactor_aij_mumps",MatGetFactor_aij_mumps);CHKERRQ(ierr); 3961 #endif 3962 #if defined(PETSC_HAVE_UMFPACK) 3963 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_umfpack_C","MatGetFactor_seqaij_umfpack",MatGetFactor_seqaij_umfpack);CHKERRQ(ierr); 3964 #endif 3965 #if defined(PETSC_HAVE_CHOLMOD) 3966 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_cholmod_C","MatGetFactor_seqaij_cholmod",MatGetFactor_seqaij_cholmod);CHKERRQ(ierr); 3967 #endif 3968 #if defined(PETSC_HAVE_LUSOL) 3969 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_lusol_C","MatGetFactor_aij_lusol",MatGetFactor_seqaij_lusol);CHKERRQ(ierr); 3970 #endif 3971 #if defined(PETSC_HAVE_CLIQUE) 3972 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_clique_C","MatGetFactor_aij_clique",MatGetFactor_aij_clique);CHKERRQ(ierr); 3973 #endif 3974 3975 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_petsc_C","MatGetFactor_seqaij_petsc",MatGetFactor_seqaij_petsc);CHKERRQ(ierr); 3976 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactorAvailable_petsc_C","MatGetFactorAvailable_seqaij_petsc",MatGetFactorAvailable_seqaij_petsc);CHKERRQ(ierr); 3977 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatGetFactor_bas_C","MatGetFactor_seqaij_bas",MatGetFactor_seqaij_bas);CHKERRQ(ierr); 3978 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetColumnIndices_C","MatSeqAIJSetColumnIndices_SeqAIJ",MatSeqAIJSetColumnIndices_SeqAIJ);CHKERRQ(ierr); 3979 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatStoreValues_C","MatStoreValues_SeqAIJ",MatStoreValues_SeqAIJ);CHKERRQ(ierr); 3980 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatRetrieveValues_C","MatRetrieveValues_SeqAIJ",MatRetrieveValues_SeqAIJ);CHKERRQ(ierr); 3981 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqsbaij_C","MatConvert_SeqAIJ_SeqSBAIJ",MatConvert_SeqAIJ_SeqSBAIJ);CHKERRQ(ierr); 3982 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqbaij_C","MatConvert_SeqAIJ_SeqBAIJ",MatConvert_SeqAIJ_SeqBAIJ);CHKERRQ(ierr); 3983 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijperm_C","MatConvert_SeqAIJ_SeqAIJPERM",MatConvert_SeqAIJ_SeqAIJPERM);CHKERRQ(ierr); 3984 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatConvert_seqaij_seqaijcrl_C","MatConvert_SeqAIJ_SeqAIJCRL",MatConvert_SeqAIJ_SeqAIJCRL);CHKERRQ(ierr); 3985 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsTranspose_C","MatIsTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3986 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatIsHermitianTranspose_C","MatIsHermitianTranspose_SeqAIJ",MatIsTranspose_SeqAIJ);CHKERRQ(ierr); 3987 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocation_C","MatSeqAIJSetPreallocation_SeqAIJ",MatSeqAIJSetPreallocation_SeqAIJ);CHKERRQ(ierr); 3988 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatSeqAIJSetPreallocationCSR_C","MatSeqAIJSetPreallocationCSR_SeqAIJ",MatSeqAIJSetPreallocationCSR_SeqAIJ);CHKERRQ(ierr); 3989 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatReorderForNonzeroDiagonal_C","MatReorderForNonzeroDiagonal_SeqAIJ",MatReorderForNonzeroDiagonal_SeqAIJ);CHKERRQ(ierr); 3990 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMult_seqdense_seqaij_C","MatMatMult_SeqDense_SeqAIJ",MatMatMult_SeqDense_SeqAIJ);CHKERRQ(ierr); 3991 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultSymbolic_seqdense_seqaij_C","MatMatMultSymbolic_SeqDense_SeqAIJ",MatMatMultSymbolic_SeqDense_SeqAIJ);CHKERRQ(ierr); 3992 ierr = PetscObjectComposeFunctionDynamic((PetscObject)B,"MatMatMultNumeric_seqdense_seqaij_C","MatMatMultNumeric_SeqDense_SeqAIJ",MatMatMultNumeric_SeqDense_SeqAIJ);CHKERRQ(ierr); 3993 ierr = MatCreate_SeqAIJ_Inode(B);CHKERRQ(ierr); 3994 ierr = PetscObjectChangeTypeName((PetscObject)B,MATSEQAIJ);CHKERRQ(ierr); 3995 PetscFunctionReturn(0); 3996 } 3997 EXTERN_C_END 3998 3999 #undef __FUNCT__ 4000 #define __FUNCT__ "MatDuplicateNoCreate_SeqAIJ" 4001 /* 4002 Given a matrix generated with MatGetFactor() duplicates all the information in A into B 4003 */ 4004 PetscErrorCode MatDuplicateNoCreate_SeqAIJ(Mat C,Mat A,MatDuplicateOption cpvalues,PetscBool mallocmatspace) 4005 { 4006 Mat_SeqAIJ *c,*a = (Mat_SeqAIJ*)A->data; 4007 PetscErrorCode ierr; 4008 PetscInt i,m = A->rmap->n; 4009 4010 PetscFunctionBegin; 4011 c = (Mat_SeqAIJ*)C->data; 4012 4013 C->factortype = A->factortype; 4014 c->row = 0; 4015 c->col = 0; 4016 c->icol = 0; 4017 c->reallocs = 0; 4018 4019 C->assembled = PETSC_TRUE; 4020 4021 ierr = PetscLayoutReference(A->rmap,&C->rmap);CHKERRQ(ierr); 4022 ierr = PetscLayoutReference(A->cmap,&C->cmap);CHKERRQ(ierr); 4023 4024 ierr = PetscMalloc2(m,PetscInt,&c->imax,m,PetscInt,&c->ilen);CHKERRQ(ierr); 4025 ierr = PetscLogObjectMemory(C, 2*m*sizeof(PetscInt));CHKERRQ(ierr); 4026 for (i=0; i<m; i++) { 4027 c->imax[i] = a->imax[i]; 4028 c->ilen[i] = a->ilen[i]; 4029 } 4030 4031 /* allocate the matrix space */ 4032 if (mallocmatspace){ 4033 ierr = PetscMalloc3(a->i[m],PetscScalar,&c->a,a->i[m],PetscInt,&c->j,m+1,PetscInt,&c->i);CHKERRQ(ierr); 4034 ierr = PetscLogObjectMemory(C, a->i[m]*(sizeof(PetscScalar)+sizeof(PetscInt))+(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4035 c->singlemalloc = PETSC_TRUE; 4036 ierr = PetscMemcpy(c->i,a->i,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4037 if (m > 0) { 4038 ierr = PetscMemcpy(c->j,a->j,(a->i[m])*sizeof(PetscInt));CHKERRQ(ierr); 4039 if (cpvalues == MAT_COPY_VALUES) { 4040 ierr = PetscMemcpy(c->a,a->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4041 } else { 4042 ierr = PetscMemzero(c->a,(a->i[m])*sizeof(PetscScalar));CHKERRQ(ierr); 4043 } 4044 } 4045 } 4046 4047 c->ignorezeroentries = a->ignorezeroentries; 4048 c->roworiented = a->roworiented; 4049 c->nonew = a->nonew; 4050 if (a->diag) { 4051 ierr = PetscMalloc((m+1)*sizeof(PetscInt),&c->diag);CHKERRQ(ierr); 4052 ierr = PetscLogObjectMemory(C,(m+1)*sizeof(PetscInt));CHKERRQ(ierr); 4053 for (i=0; i<m; i++) { 4054 c->diag[i] = a->diag[i]; 4055 } 4056 } else c->diag = 0; 4057 c->solve_work = 0; 4058 c->saved_values = 0; 4059 c->idiag = 0; 4060 c->ssor_work = 0; 4061 c->keepnonzeropattern = a->keepnonzeropattern; 4062 c->free_a = PETSC_TRUE; 4063 c->free_ij = PETSC_TRUE; 4064 c->xtoy = 0; 4065 c->XtoY = 0; 4066 4067 c->rmax = a->rmax; 4068 c->nz = a->nz; 4069 c->maxnz = a->nz; /* Since we allocate exactly the right amount */ 4070 C->preallocated = PETSC_TRUE; 4071 4072 c->compressedrow.use = a->compressedrow.use; 4073 c->compressedrow.nrows = a->compressedrow.nrows; 4074 c->compressedrow.check = a->compressedrow.check; 4075 if (a->compressedrow.use){ 4076 i = a->compressedrow.nrows; 4077 ierr = PetscMalloc2(i+1,PetscInt,&c->compressedrow.i,i,PetscInt,&c->compressedrow.rindex);CHKERRQ(ierr); 4078 ierr = PetscMemcpy(c->compressedrow.i,a->compressedrow.i,(i+1)*sizeof(PetscInt));CHKERRQ(ierr); 4079 ierr = PetscMemcpy(c->compressedrow.rindex,a->compressedrow.rindex,i*sizeof(PetscInt));CHKERRQ(ierr); 4080 } else { 4081 c->compressedrow.use = PETSC_FALSE; 4082 c->compressedrow.i = PETSC_NULL; 4083 c->compressedrow.rindex = PETSC_NULL; 4084 } 4085 C->same_nonzero = A->same_nonzero; 4086 ierr = MatDuplicate_SeqAIJ_Inode(A,cpvalues,&C);CHKERRQ(ierr); 4087 4088 ierr = PetscFListDuplicate(((PetscObject)A)->qlist,&((PetscObject)C)->qlist);CHKERRQ(ierr); 4089 PetscFunctionReturn(0); 4090 } 4091 4092 #undef __FUNCT__ 4093 #define __FUNCT__ "MatDuplicate_SeqAIJ" 4094 PetscErrorCode MatDuplicate_SeqAIJ(Mat A,MatDuplicateOption cpvalues,Mat *B) 4095 { 4096 PetscErrorCode ierr; 4097 4098 PetscFunctionBegin; 4099 ierr = MatCreate(((PetscObject)A)->comm,B);CHKERRQ(ierr); 4100 ierr = MatSetSizes(*B,A->rmap->n,A->cmap->n,A->rmap->n,A->cmap->n);CHKERRQ(ierr); 4101 ierr = MatSetBlockSizes(*B,A->rmap->bs,A->cmap->bs);CHKERRQ(ierr); 4102 ierr = MatSetType(*B,((PetscObject)A)->type_name);CHKERRQ(ierr); 4103 ierr = MatDuplicateNoCreate_SeqAIJ(*B,A,cpvalues,PETSC_TRUE);CHKERRQ(ierr); 4104 PetscFunctionReturn(0); 4105 } 4106 4107 #undef __FUNCT__ 4108 #define __FUNCT__ "MatLoad_SeqAIJ" 4109 PetscErrorCode MatLoad_SeqAIJ(Mat newMat, PetscViewer viewer) 4110 { 4111 Mat_SeqAIJ *a; 4112 PetscErrorCode ierr; 4113 PetscInt i,sum,nz,header[4],*rowlengths = 0,M,N,rows,cols; 4114 int fd; 4115 PetscMPIInt size; 4116 MPI_Comm comm; 4117 PetscInt bs = 1; 4118 4119 PetscFunctionBegin; 4120 ierr = PetscObjectGetComm((PetscObject)viewer,&comm);CHKERRQ(ierr); 4121 ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr); 4122 if (size > 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"view must have one processor"); 4123 4124 ierr = PetscOptionsBegin(comm,PETSC_NULL,"Options for loading SEQAIJ matrix","Mat");CHKERRQ(ierr); 4125 ierr = PetscOptionsInt("-matload_block_size","Set the blocksize used to store the matrix","MatLoad",bs,&bs,PETSC_NULL);CHKERRQ(ierr); 4126 ierr = PetscOptionsEnd();CHKERRQ(ierr); 4127 4128 ierr = PetscViewerBinaryGetDescriptor(viewer,&fd);CHKERRQ(ierr); 4129 ierr = PetscBinaryRead(fd,header,4,PETSC_INT);CHKERRQ(ierr); 4130 if (header[0] != MAT_FILE_CLASSID) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"not matrix object in file"); 4131 M = header[1]; N = header[2]; nz = header[3]; 4132 4133 if (nz < 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED,"Matrix stored in special format on disk,cannot load as SeqAIJ"); 4134 4135 /* read in row lengths */ 4136 ierr = PetscMalloc(M*sizeof(PetscInt),&rowlengths);CHKERRQ(ierr); 4137 ierr = PetscBinaryRead(fd,rowlengths,M,PETSC_INT);CHKERRQ(ierr); 4138 4139 /* check if sum of rowlengths is same as nz */ 4140 for (i=0,sum=0; i< M; i++) sum +=rowlengths[i]; 4141 if (sum != nz) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Inconsistant matrix data in file. no-nonzeros = %d, sum-row-lengths = %d\n",nz,sum); 4142 4143 /* set global size if not set already*/ 4144 if (newMat->rmap->n < 0 && newMat->rmap->N < 0 && newMat->cmap->n < 0 && newMat->cmap->N < 0) { 4145 ierr = MatSetSizes(newMat,PETSC_DECIDE,PETSC_DECIDE,M,N);CHKERRQ(ierr); 4146 } else { 4147 /* if sizes and type are already set, check if the vector global sizes are correct */ 4148 ierr = MatGetSize(newMat,&rows,&cols);CHKERRQ(ierr); 4149 if (rows < 0 && cols < 0){ /* user might provide local size instead of global size */ 4150 ierr = MatGetLocalSize(newMat,&rows,&cols);CHKERRQ(ierr); 4151 } 4152 if (M != rows || N != cols) SETERRQ4(PETSC_COMM_SELF,PETSC_ERR_FILE_UNEXPECTED, "Matrix in file of different length (%d, %d) than the input matrix (%d, %d)",M,N,rows,cols); 4153 } 4154 ierr = MatSeqAIJSetPreallocation_SeqAIJ(newMat,0,rowlengths);CHKERRQ(ierr); 4155 a = (Mat_SeqAIJ*)newMat->data; 4156 4157 ierr = PetscBinaryRead(fd,a->j,nz,PETSC_INT);CHKERRQ(ierr); 4158 4159 /* read in nonzero values */ 4160 ierr = PetscBinaryRead(fd,a->a,nz,PETSC_SCALAR);CHKERRQ(ierr); 4161 4162 /* set matrix "i" values */ 4163 a->i[0] = 0; 4164 for (i=1; i<= M; i++) { 4165 a->i[i] = a->i[i-1] + rowlengths[i-1]; 4166 a->ilen[i-1] = rowlengths[i-1]; 4167 } 4168 ierr = PetscFree(rowlengths);CHKERRQ(ierr); 4169 4170 ierr = MatAssemblyBegin(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4171 ierr = MatAssemblyEnd(newMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4172 if (bs > 1) {ierr = MatSetBlockSize(newMat,bs);CHKERRQ(ierr);} 4173 PetscFunctionReturn(0); 4174 } 4175 4176 #undef __FUNCT__ 4177 #define __FUNCT__ "MatEqual_SeqAIJ" 4178 PetscErrorCode MatEqual_SeqAIJ(Mat A,Mat B,PetscBool * flg) 4179 { 4180 Mat_SeqAIJ *a = (Mat_SeqAIJ *)A->data,*b = (Mat_SeqAIJ *)B->data; 4181 PetscErrorCode ierr; 4182 #if defined(PETSC_USE_COMPLEX) 4183 PetscInt k; 4184 #endif 4185 4186 PetscFunctionBegin; 4187 /* If the matrix dimensions are not equal,or no of nonzeros */ 4188 if ((A->rmap->n != B->rmap->n) || (A->cmap->n != B->cmap->n) ||(a->nz != b->nz)) { 4189 *flg = PETSC_FALSE; 4190 PetscFunctionReturn(0); 4191 } 4192 4193 /* if the a->i are the same */ 4194 ierr = PetscMemcmp(a->i,b->i,(A->rmap->n+1)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4195 if (!*flg) PetscFunctionReturn(0); 4196 4197 /* if a->j are the same */ 4198 ierr = PetscMemcmp(a->j,b->j,(a->nz)*sizeof(PetscInt),flg);CHKERRQ(ierr); 4199 if (!*flg) PetscFunctionReturn(0); 4200 4201 /* if a->a are the same */ 4202 #if defined(PETSC_USE_COMPLEX) 4203 for (k=0; k<a->nz; k++){ 4204 if (PetscRealPart(a->a[k]) != PetscRealPart(b->a[k]) || PetscImaginaryPart(a->a[k]) != PetscImaginaryPart(b->a[k])){ 4205 *flg = PETSC_FALSE; 4206 PetscFunctionReturn(0); 4207 } 4208 } 4209 #else 4210 ierr = PetscMemcmp(a->a,b->a,(a->nz)*sizeof(PetscScalar),flg);CHKERRQ(ierr); 4211 #endif 4212 PetscFunctionReturn(0); 4213 } 4214 4215 #undef __FUNCT__ 4216 #define __FUNCT__ "MatCreateSeqAIJWithArrays" 4217 /*@ 4218 MatCreateSeqAIJWithArrays - Creates an sequential AIJ matrix using matrix elements (in CSR format) 4219 provided by the user. 4220 4221 Collective on MPI_Comm 4222 4223 Input Parameters: 4224 + comm - must be an MPI communicator of size 1 4225 . m - number of rows 4226 . n - number of columns 4227 . i - row indices 4228 . j - column indices 4229 - a - matrix values 4230 4231 Output Parameter: 4232 . mat - the matrix 4233 4234 Level: intermediate 4235 4236 Notes: 4237 The i, j, and a arrays are not copied by this routine, the user must free these arrays 4238 once the matrix is destroyed and not before 4239 4240 You cannot set new nonzero locations into this matrix, that will generate an error. 4241 4242 The i and j indices are 0 based 4243 4244 The format which is used for the sparse matrix input, is equivalent to a 4245 row-major ordering.. i.e for the following matrix, the input data expected is 4246 as shown: 4247 4248 1 0 0 4249 2 0 3 4250 4 5 6 4251 4252 i = {0,1,3,6} [size = nrow+1 = 3+1] 4253 j = {0,0,2,0,1,2} [size = nz = 6]; values must be sorted for each row 4254 v = {1,2,3,4,5,6} [size = nz = 6] 4255 4256 4257 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateMPIAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4258 4259 @*/ 4260 PetscErrorCode MatCreateSeqAIJWithArrays(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat) 4261 { 4262 PetscErrorCode ierr; 4263 PetscInt ii; 4264 Mat_SeqAIJ *aij; 4265 #if defined(PETSC_USE_DEBUG) 4266 PetscInt jj; 4267 #endif 4268 4269 PetscFunctionBegin; 4270 if (i[0]) { 4271 SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"i (row indices) must start with 0"); 4272 } 4273 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4274 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4275 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4276 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4277 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,MAT_SKIP_ALLOCATION,0);CHKERRQ(ierr); 4278 aij = (Mat_SeqAIJ*)(*mat)->data; 4279 ierr = PetscMalloc2(m,PetscInt,&aij->imax,m,PetscInt,&aij->ilen);CHKERRQ(ierr); 4280 4281 aij->i = i; 4282 aij->j = j; 4283 aij->a = a; 4284 aij->singlemalloc = PETSC_FALSE; 4285 aij->nonew = -1; /*this indicates that inserting a new value in the matrix that generates a new nonzero is an error*/ 4286 aij->free_a = PETSC_FALSE; 4287 aij->free_ij = PETSC_FALSE; 4288 4289 for (ii=0; ii<m; ii++) { 4290 aij->ilen[ii] = aij->imax[ii] = i[ii+1] - i[ii]; 4291 #if defined(PETSC_USE_DEBUG) 4292 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]); 4293 for (jj=i[ii]+1; jj<i[ii+1]; jj++) { 4294 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); 4295 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); 4296 } 4297 #endif 4298 } 4299 #if defined(PETSC_USE_DEBUG) 4300 for (ii=0; ii<aij->i[m]; ii++) { 4301 if (j[ii] < 0) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Negative column index at location = %d index = %d",ii,j[ii]); 4302 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]); 4303 } 4304 #endif 4305 4306 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4307 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4308 PetscFunctionReturn(0); 4309 } 4310 #undef __FUNCT__ 4311 #define __FUNCT__ "MatCreateSeqAIJFromTriple" 4312 /*@C 4313 MatCreateSeqAIJFromTriple - Creates an sequential AIJ matrix using matrix elements (in COO format) 4314 provided by the user. 4315 4316 Collective on MPI_Comm 4317 4318 Input Parameters: 4319 + comm - must be an MPI communicator of size 1 4320 . m - number of rows 4321 . n - number of columns 4322 . i - row indices 4323 . j - column indices 4324 . a - matrix values 4325 . nz - number of nonzeros 4326 - idx - 0 or 1 based 4327 4328 Output Parameter: 4329 . mat - the matrix 4330 4331 Level: intermediate 4332 4333 Notes: 4334 The i and j indices are 0 based 4335 4336 The format which is used for the sparse matrix input, is equivalent to a 4337 row-major ordering.. i.e for the following matrix, the input data expected is 4338 as shown: 4339 4340 1 0 0 4341 2 0 3 4342 4 5 6 4343 4344 i = {0,1,1,2,2,2} 4345 j = {0,0,2,0,1,2} 4346 v = {1,2,3,4,5,6} 4347 4348 4349 .seealso: MatCreate(), MatCreateAIJ(), MatCreateSeqAIJ(), MatCreateSeqAIJWithArrays(), MatMPIAIJSetPreallocationCSR() 4350 4351 @*/ 4352 PetscErrorCode MatCreateSeqAIJFromTriple(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt* i,PetscInt*j,PetscScalar *a,Mat *mat,PetscInt nz,PetscBool idx) 4353 { 4354 PetscErrorCode ierr; 4355 PetscInt ii, *nnz, one = 1,row,col; 4356 4357 4358 PetscFunctionBegin; 4359 ierr = PetscMalloc(m*sizeof(PetscInt),&nnz);CHKERRQ(ierr); 4360 ierr = PetscMemzero(nnz,m*sizeof(PetscInt));CHKERRQ(ierr); 4361 for (ii = 0; ii < nz; ii++){ 4362 nnz[i[ii]] += 1; 4363 } 4364 ierr = MatCreate(comm,mat);CHKERRQ(ierr); 4365 ierr = MatSetSizes(*mat,m,n,m,n);CHKERRQ(ierr); 4366 /* ierr = MatSetBlockSizes(*mat,,);CHKERRQ(ierr); */ 4367 ierr = MatSetType(*mat,MATSEQAIJ);CHKERRQ(ierr); 4368 ierr = MatSeqAIJSetPreallocation_SeqAIJ(*mat,0,nnz);CHKERRQ(ierr); 4369 for (ii = 0; ii < nz; ii++){ 4370 if (idx){ 4371 row = i[ii] - 1; 4372 col = j[ii] - 1; 4373 } else { 4374 row = i[ii]; 4375 col = j[ii]; 4376 } 4377 ierr = MatSetValues(*mat,one,&row,one,&col,&a[ii],ADD_VALUES);CHKERRQ(ierr); 4378 } 4379 ierr = MatAssemblyBegin(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4380 ierr = MatAssemblyEnd(*mat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 4381 ierr = PetscFree(nnz);CHKERRQ(ierr); 4382 PetscFunctionReturn(0); 4383 } 4384 4385 #undef __FUNCT__ 4386 #define __FUNCT__ "MatSetColoring_SeqAIJ" 4387 PetscErrorCode MatSetColoring_SeqAIJ(Mat A,ISColoring coloring) 4388 { 4389 PetscErrorCode ierr; 4390 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4391 4392 PetscFunctionBegin; 4393 if (coloring->ctype == IS_COLORING_GLOBAL) { 4394 ierr = ISColoringReference(coloring);CHKERRQ(ierr); 4395 a->coloring = coloring; 4396 } else if (coloring->ctype == IS_COLORING_GHOSTED) { 4397 PetscInt i,*larray; 4398 ISColoring ocoloring; 4399 ISColoringValue *colors; 4400 4401 /* set coloring for diagonal portion */ 4402 ierr = PetscMalloc(A->cmap->n*sizeof(PetscInt),&larray);CHKERRQ(ierr); 4403 for (i=0; i<A->cmap->n; i++) { 4404 larray[i] = i; 4405 } 4406 ierr = ISGlobalToLocalMappingApply(A->cmap->mapping,IS_GTOLM_MASK,A->cmap->n,larray,PETSC_NULL,larray);CHKERRQ(ierr); 4407 ierr = PetscMalloc(A->cmap->n*sizeof(ISColoringValue),&colors);CHKERRQ(ierr); 4408 for (i=0; i<A->cmap->n; i++) { 4409 colors[i] = coloring->colors[larray[i]]; 4410 } 4411 ierr = PetscFree(larray);CHKERRQ(ierr); 4412 ierr = ISColoringCreate(PETSC_COMM_SELF,coloring->n,A->cmap->n,colors,&ocoloring);CHKERRQ(ierr); 4413 a->coloring = ocoloring; 4414 } 4415 PetscFunctionReturn(0); 4416 } 4417 4418 #undef __FUNCT__ 4419 #define __FUNCT__ "MatSetValuesAdifor_SeqAIJ" 4420 PetscErrorCode MatSetValuesAdifor_SeqAIJ(Mat A,PetscInt nl,void *advalues) 4421 { 4422 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4423 PetscInt m = A->rmap->n,*ii = a->i,*jj = a->j,nz,i,j; 4424 MatScalar *v = a->a; 4425 PetscScalar *values = (PetscScalar *)advalues; 4426 ISColoringValue *color; 4427 4428 PetscFunctionBegin; 4429 if (!a->coloring) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"Coloring not set for matrix"); 4430 color = a->coloring->colors; 4431 /* loop over rows */ 4432 for (i=0; i<m; i++) { 4433 nz = ii[i+1] - ii[i]; 4434 /* loop over columns putting computed value into matrix */ 4435 for (j=0; j<nz; j++) { 4436 *v++ = values[color[*jj++]]; 4437 } 4438 values += nl; /* jump to next row of derivatives */ 4439 } 4440 PetscFunctionReturn(0); 4441 } 4442 4443 /* 4444 Special version for direct calls from Fortran 4445 */ 4446 #include <petsc-private/fortranimpl.h> 4447 #if defined(PETSC_HAVE_FORTRAN_CAPS) 4448 #define matsetvaluesseqaij_ MATSETVALUESSEQAIJ 4449 #elif !defined(PETSC_HAVE_FORTRAN_UNDERSCORE) 4450 #define matsetvaluesseqaij_ matsetvaluesseqaij 4451 #endif 4452 4453 /* Change these macros so can be used in void function */ 4454 #undef CHKERRQ 4455 #define CHKERRQ(ierr) CHKERRABORT(((PetscObject)A)->comm,ierr) 4456 #undef SETERRQ2 4457 #define SETERRQ2(comm,ierr,b,c,d) CHKERRABORT(comm,ierr) 4458 #undef SETERRQ3 4459 #define SETERRQ3(comm,ierr,b,c,d,e) CHKERRABORT(comm,ierr) 4460 4461 EXTERN_C_BEGIN 4462 #undef __FUNCT__ 4463 #define __FUNCT__ "matsetvaluesseqaij_" 4464 void PETSC_STDCALL matsetvaluesseqaij_(Mat *AA,PetscInt *mm,const PetscInt im[],PetscInt *nn,const PetscInt in[],const PetscScalar v[],InsertMode *isis, PetscErrorCode *_ierr) 4465 { 4466 Mat A = *AA; 4467 PetscInt m = *mm, n = *nn; 4468 InsertMode is = *isis; 4469 Mat_SeqAIJ *a = (Mat_SeqAIJ*)A->data; 4470 PetscInt *rp,k,low,high,t,ii,row,nrow,i,col,l,rmax,N; 4471 PetscInt *imax,*ai,*ailen; 4472 PetscErrorCode ierr; 4473 PetscInt *aj,nonew = a->nonew,lastcol = -1; 4474 MatScalar *ap,value,*aa; 4475 PetscBool ignorezeroentries = a->ignorezeroentries; 4476 PetscBool roworiented = a->roworiented; 4477 4478 PetscFunctionBegin; 4479 MatCheckPreallocated(A,1); 4480 imax = a->imax; 4481 ai = a->i; 4482 ailen = a->ilen; 4483 aj = a->j; 4484 aa = a->a; 4485 4486 for (k=0; k<m; k++) { /* loop over added rows */ 4487 row = im[k]; 4488 if (row < 0) continue; 4489 #if defined(PETSC_USE_DEBUG) 4490 if (row >= A->rmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Row too large"); 4491 #endif 4492 rp = aj + ai[row]; ap = aa + ai[row]; 4493 rmax = imax[row]; nrow = ailen[row]; 4494 low = 0; 4495 high = nrow; 4496 for (l=0; l<n; l++) { /* loop over added columns */ 4497 if (in[l] < 0) continue; 4498 #if defined(PETSC_USE_DEBUG) 4499 if (in[l] >= A->cmap->n) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Column too large"); 4500 #endif 4501 col = in[l]; 4502 if (roworiented) { 4503 value = v[l + k*n]; 4504 } else { 4505 value = v[k + l*m]; 4506 } 4507 if (value == 0.0 && ignorezeroentries && (is == ADD_VALUES)) continue; 4508 4509 if (col <= lastcol) low = 0; else high = nrow; 4510 lastcol = col; 4511 while (high-low > 5) { 4512 t = (low+high)/2; 4513 if (rp[t] > col) high = t; 4514 else low = t; 4515 } 4516 for (i=low; i<high; i++) { 4517 if (rp[i] > col) break; 4518 if (rp[i] == col) { 4519 if (is == ADD_VALUES) ap[i] += value; 4520 else ap[i] = value; 4521 goto noinsert; 4522 } 4523 } 4524 if (value == 0.0 && ignorezeroentries) goto noinsert; 4525 if (nonew == 1) goto noinsert; 4526 if (nonew == -1) SETERRABORT(((PetscObject)A)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Inserting a new nonzero in the matrix"); 4527 MatSeqXAIJReallocateAIJ(A,A->rmap->n,1,nrow,row,col,rmax,aa,ai,aj,rp,ap,imax,nonew,MatScalar); 4528 N = nrow++ - 1; a->nz++; high++; 4529 /* shift up all the later entries in this row */ 4530 for (ii=N; ii>=i; ii--) { 4531 rp[ii+1] = rp[ii]; 4532 ap[ii+1] = ap[ii]; 4533 } 4534 rp[i] = col; 4535 ap[i] = value; 4536 noinsert:; 4537 low = i + 1; 4538 } 4539 ailen[row] = nrow; 4540 } 4541 A->same_nonzero = PETSC_FALSE; 4542 PetscFunctionReturnVoid(); 4543 } 4544 EXTERN_C_END 4545 4546