1. petsc
  2. PETSc
  3. petsc

Commits

Shrirang Abhyankar  committed 142808b

[petsc-maint #106790] Fixed bug in MatMult_SeqSBAIJ_1 versions for symmetric matrices having empty diagonal.

Hg-commit: 3b022781a8df6d13f2a1e23f284add510c79dace

  • Participants
  • Parent commits 17803ae
  • Branches master

Comments (0)

Files changed (1)

File src/mat/impls/sbaij/seq/relax.h

View file
     if (!nz) continue; /* Move to the next row if the current row is empty */
     nonzerorow++;
     x1   = x[i];
-    sum  = v[0]*x1;          /* diagonal term */
-    for (j=1; j<nz; j++) {
+    sum = 0;
+    for (j=0; j<nz; j++) {
       ibt  = ib[j];
       vj   = v[j];
-      sum += vj * x[ibt];   /* (strict upper triangular part of A)*x  */
-      z[ibt] += PetscConj(v[j]) * x1;    /* (strict lower triangular part of A)*x  */
+      sum += vj * x[ibt];   /* (Upper triangular part of A)*x  */
+      if(ibt != i) z[ibt] += PetscConj(v[j]) * x1;    /* (strict lower triangular part of A)*x  */
     }
     z[i] += sum;
     v    += nz;
     if (!nz) continue; /* Move to the next row if the current row is empty */
     nonzerorow++;
     x1   = x[i];
-    sum  = v[0]*x1;                /* diagonal term */
+    sum = 0;
     PetscPrefetchBlock(ib+nz,nz,0,PETSC_PREFETCH_HINT_NTA); /* Indices for the next row (assumes same size as this one) */
     PetscPrefetchBlock(v+nz,nz,0,PETSC_PREFETCH_HINT_NTA);  /* Entries for the next row */
-    for (j=1; j<nz; j++) {
+    for (j=0; j<nz; j++) {
       ibt = ib[j];
       vj  = v[j];
-      z[ibt] += vj * x1;       /* (strict lower triangular part of A)*x  */
-      sum    += vj * x[ibt]; /* (strict upper triangular part of A)*x  */
+      if(ibt != i) z[ibt] += vj * x1;       /* (strict lower triangular part of A)*x  */
+      sum    += vj * x[ibt]; /* (Upper triangular part of A)*x  */
     }
     z[i] += sum;
     v    += nz;