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