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