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