SCM

SCM Repository

[matrix] Diff of /pkg/src/dgBCMatrix.c
ViewVC logotype

Diff of /pkg/src/dgBCMatrix.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 552, Tue Feb 15 23:39:40 2005 UTC revision 553, Wed Feb 16 13:43:17 2005 UTC
# Line 180  Line 180 
180    
181                  if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, ii);                  if (K < 0) error("cscb_syrk: C[%d,%d] not defined", ii, ii);
182                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];                  if (scalar) Cx[K] += alpha * Ax[k] * Ax[k];
183                  else F77_CALL(dsyrk)((uplo == UPP)?"U":"L", "N", cdims, adims + 1,                  else F77_CALL(dsyrk)((uplo == UPP) ? "U" : "L", "N",
184                                         cdims, adims + 1,
185                                       &alpha, Ax + k * asz, adims,                                       &alpha, Ax + k * asz, adims,
186                                       &one, Cx + K * csz, cdims);                                       &one, Cx + K * csz, cdims);
187                  for (kk = k+1; kk < k2; kk++) {                  for (kk = k+1; kk < k2; kk++) {
# Line 189  Line 190 
190                          check_csc_index(Cp, Ci, jj, ii, 0);                          check_csc_index(Cp, Ci, jj, ii, 0);
191    
192                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];                      if (scalar) Cx[K] += alpha * Ax[k] * Ax[kk];
193                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1, adims + 1,                      else F77_CALL(dgemm)("N", "T", cdims, cdims + 1,
194                                           &alpha, Ax + ((uplo==UPP)?k:kk) * asz, adims,                                           adims + 1, &alpha,
195                                             Ax+((uplo==UPP)?k:kk)*asz, adims,
196                                           Ax + ((uplo==UPP)?kk:k) * asz, adims,                                           Ax + ((uplo==UPP)?kk:k) * asz, adims,
197                                           &one, Cx + K * csz, cdims);                                           &one, Cx + K * csz, cdims);
198                  }                  }
# Line 269  Line 271 
271              p2 = Ap[k+1];              p2 = Ap[k+1];
272              for (p = Ap[k]; p < p2; p++) {              for (p = Ap[k]; p < p2; p++) {
273                  i = Ai[p];      /* get A[i,k] */                  i = Ai[p];      /* get A[i,k] */
274                  if (i <= k) {   /* [i,k] in upper triangle? Should always be true */                  if (i <= k) {   /* [i,k] in upper triangle? Should be true */
275                                  /* copy A(i,k) into Y */                                  /* copy A(i,k) into Y */
276                      Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);                      Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);
277                      /* follow path from i to root of etree, stop at flagged node */                      /* follow path from i to root of etree,
278                         * stop at flagged node */
279                      for (len = 0; Flag[i] != k; i = Parent[i]) {                      for (len = 0; Flag[i] != k; i = Parent[i]) {
280                          Pattern[len++] = i; /* L[k,i] is nonzero */                          Pattern[len++] = i; /* L[k,i] is nonzero */
281                          Flag[i] = k; /* mark i as visited */                          Flag[i] = k; /* mark i as visited */
# Line 283  Line 286 
286                  }                  }
287              }              }
288              /* Pattern [top ... n-1] now contains nonzero pattern of L[,k] */              /* Pattern [top ... n-1] now contains nonzero pattern of L[,k] */
289              /* compute numerical values in kth row of L (a sparse triangular solve) */              /* compute numerical values in kth row of L
290                 * (a sparse triangular solve) */
291              Memcpy(Dx + k * ncisqr, Y + k * ncisqr, ncisqr); /* get D[,,k] */              Memcpy(Dx + k * ncisqr, Y + k * ncisqr, ncisqr); /* get D[,,k] */
292              AZERO(Y + k*ncisqr, ncisqr); /* clear Y[,,k] */              AZERO(Y + k*ncisqr, ncisqr); /* clear Y[,,k] */
293              for (; top < n; top++) {              for (; top < n; top++) {
# Line 296  Line 300 
300                                      Lx + p*ncisqr, &nci, Yi, &nci,                                      Lx + p*ncisqr, &nci, Yi, &nci,
301                                      &one, Y + Li[p]*ncisqr, &nci);                                      &one, Y + Li[p]*ncisqr, &nci);
302                  }                  }
303                  /* FIXME: Check that this is the correct order and transposition */                  /* FIXME: Is this the correct order and transposition? */
304                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,
305                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
306                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,

Legend:
Removed from v.552  
changed lines
  Added in v.553

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge