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