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