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