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