# SCM Repository

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

# Diff of /pkg/src/dgBCMatrix.c

revision 553, Wed Feb 16 13:43:17 2005 UTC revision 556, Wed Feb 16 20:49:44 2005 UTC
# Line 1  Line 1
1  #include "dgBCMatrix.h"  #include "dgBCMatrix.h"
2  /* TODO  /* TODO
3   *  - code for trans = 'T' in cscb_syrk   *  - code for trans = 'T' in cscb_syrk
4   *  - code for non-trivial cscb_trmm and cscb_ldl   *  - code for non-trivial cscb_trmm
5   */   */
6
7  SEXP dgBCMatrix_validate(SEXP x)  SEXP dgBCMatrix_validate(SEXP x)
# Line 201  Line 201
201      }      }
202  }  }
203
204    static void
205    copy_transpose(double dest[], const double src[], int n)
206    {
207        int i, j;
208        for (i = 0; i < n; i++) {
209            for (j = 0; j < n; j++) {
210                dest[i + j * n] = src[j + i * n];
211            }
212        }
213    }
214
215  /**  /**
216   * Create the LD^{T/2}D^{1/2}L' decomposition of the positive definite   * Create the LD^{T/2}D^{1/2}L' decomposition of the positive definite
217   * symmetric dgBCMatrix matrix A (upper triangle stored) in L and D^{1/2}.   * symmetric dgBCMatrix matrix A (upper triangle stored) in L and D^{1/2}.
# Line 241  Line 252
252          Memcpy(Dx, Ax, ncisqr * n);          Memcpy(Dx, Ax, ncisqr * n);
253          for (j = 0; j < n; j++) { /* form D_i^{1/2} */          for (j = 0; j < n; j++) { /* form D_i^{1/2} */
254              F77_CALL(dpotrf)("U", &nci, Dx + j * ncisqr, &nci, &k);              F77_CALL(dpotrf)("U", &nci, Dx + j * ncisqr, &nci, &k);
255              /* Should we just return k instead of signaling an error? */              if (k) return j; /* return block number, not col no. */
/* The caller should check for return value < n */
if (k) error("D[ , , %d], level %d, is not positive definite",
j + 1);
256          }          }
257          return n;          return n;
258      }      }
259      if (nci == 1) {      if (nci == 1) {
260          k = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent, Li, Lx, Dx,          k = R_ldl_numeric(n, Ap, Ai, Ax, Lp, Parent, Li, Lx, Dx,
261                            (int *) NULL, (int *) NULL);                            (int *) NULL, (int *) NULL);
262          if (k < n) error("cscb_ldl: error code %k from R_ldl_numeric", k);          if (k < n) return k;
263          for (j = 0; j < n; j++) Dx[j] = sqrt(Dx[j]);          for (j = 0; j < n; j++) Dx[j] = sqrt(Dx[j]);
264          return n;          return n;
265      } else {               /* Copy of ldl_numeric from the LDL package      } else {               /* Copy of ldl_numeric from the LDL package
# Line 271  Line 279
279              p2 = Ap[k+1];              p2 = Ap[k+1];
280              for (p = Ap[k]; p < p2; p++) {              for (p = Ap[k]; p < p2; p++) {
281                  i = Ai[p];      /* get A[i,k] */                  i = Ai[p];      /* get A[i,k] */
282                  if (i <= k) {   /* [i,k] in upper triangle? Should be true */                  if (i > k) error("cscb_ldl: A has nonzeros below diagonal");
283                                  /* copy A(i,k) into Y */                                  /* copy A(i,k) into Y */
284                      Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);                      Memcpy(Y + i * ncisqr, Ax + p * ncisqr, ncisqr);
285                      /* follow path from i to root of etree,                      /* follow path from i to root of etree,
# Line 284  Line 292
292                          Pattern[--top] = Pattern[--len];                          Pattern[--top] = Pattern[--len];
293                      }                      }
294                  }                  }
}
295              /* Pattern [top ... n-1] now contains nonzero pattern of L[,k] */              /* Pattern [top ... n-1] now contains nonzero pattern of L[,k] */
296              /* compute numerical values in kth row of L              /* compute numerical values in kth row of L
297               * (a sparse triangular solve) */               * (a sparse triangular solve) */
# Line 301  Line 308
308                                      &one, Y + Li[p]*ncisqr, &nci);                                      &one, Y + Li[p]*ncisqr, &nci);
309                  }                  }
310                  /* FIXME: Is this the correct order and transposition? */                  /* FIXME: Is this the correct order and transposition? */
311                  F77_CALL(dtrsm)("R", "U", "T", "N", &nci, &nci,                  F77_CALL(dtrsm)("L", "U", "T", "N", &nci, &nci,
312                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
313                  F77_CALL(dsyrk)("U", "N", &nci, &nci, &minus1, Yi, &nci,                  F77_CALL(dsyrk)("U", "T", &nci, &nci, &minus1, Yi, &nci,
314                                  &one, Dx + k*ncisqr, &nci);                                  &one, Dx + k*ncisqr, &nci);
315                  F77_CALL(dtrsm)("R", "U", "N", "N", &nci, &nci,                  F77_CALL(dtrsm)("L", "U", "N", "N", &nci, &nci,
316                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);                                  &one, Dx + i*ncisqr, &nci, Yi, &nci);
317                  Li[p] = k;      /* store L[k,i] in column form of L */                  Li[p] = k;      /* store L[k,i] in column form of L */
318                  Memcpy(Lx + p * ncisqr, Yi, ncisqr);                  /* Yi contains L[k,i]', not L[k,i] */
319                    copy_transpose(Lx + p * ncisqr, Yi, nci);
320                  Lnz[i]++;       /* increment count of nonzeros in col i */                  Lnz[i]++;       /* increment count of nonzeros in col i */
321              }              }
322              F77_CALL(dpotrf)("U", &nci, Dx + k*ncisqr, &nci, &info);              F77_CALL(dpotrf)("U", &nci, Dx + k*ncisqr, &nci, &info);
323              if (info) {              if (info) {
324                  Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);                  Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);
325                  return k;           /* failure, D[,,k] is not positive definite */                  return k;  /* failure, D[,,k] not positive definite */
326              }              }
327          }          }
328          Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);          Free(Y); Free(Yi); Free(Pattern); Free(Flag); Free(Lnz);
# Line 375  Line 383
383   * @param ldb leading dimension of B as declared in the calling function   * @param ldb leading dimension of B as declared in the calling function
384   */   */
385  void  void
386  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa, enum CBLAS_DIAG diag,  cscb_trsm(enum CBLAS_UPLO uplo, enum CBLAS_TRANSPOSE transa,
387            double alpha, SEXP A, double B[], int m, int n, int ldb)            enum CBLAS_DIAG diag, double alpha, SEXP A,
388              double B[], int m, int n, int ldb)
389  {  {
390      SEXP ApP = GET_SLOT(A, Matrix_pSym),      SEXP ApP = GET_SLOT(A, Matrix_pSym),
391          AxP = GET_SLOT(A, Matrix_xSym);          AxP = GET_SLOT(A, Matrix_xSym);
# Line 542  Line 551
551              double *BTx = Calloc(nnz, double), *rhs;              double *BTx = Calloc(nnz, double), *rhs;
552
553                                  /* transpose B */                                  /* transpose B */
554              for (i = 0, nrbB = -1; i < nnz; i++) if (Bi[i] > nrbB) nrbB = Bi[i];              for (i = 0, nrbB = -1; i < nnz; i++)
555                    if (Bi[i] > nrbB) nrbB = Bi[i];
556              BTp = Calloc(nrbB, int);              BTp = Calloc(nrbB, int);
557              triplet_to_col(ncbB, nrbB, nnz, tmp, Bi, Bx, BTp, BTi, BTx);              triplet_to_col(ncbB, nrbB, nnz, tmp, Bi, Bx, BTp, BTi, BTx);
558                                  /* sanity check */                                  /* sanity check */
# Line 552  Line 562
562              rhs = Calloc(ncbB, double);              rhs = Calloc(ncbB, double);
563              AZERO(Bx, nnz);     /* zero the result */              AZERO(Bx, nnz);     /* zero the result */
564              for (i = 0; i < nrbB; i++) {              for (i = 0; i < nrbB; i++) {
565                  R_ldl_lsolve(ncbB, expand_csc_column(rhs, ncbB, i, BTp, BTi, BTx),                  R_ldl_lsolve(ncbB,
566                                 expand_csc_column(rhs, ncbB, i, BTp, BTi, BTx),
567                               Ap, Ai, Ax);                               Ap, Ai, Ax);
568                  for (j = 0; j < ncbB; j++) { /* write non-zeros in sol'n into B */                  /* write non-zeros in sol'n into B */
569                    for (j = 0; j < ncbB; j++) {
570                      if (BTx[j]) Bx[check_csc_index(Bp, Bi, j, i, 0)] = BTx[j];                      if (BTx[j]) Bx[check_csc_index(Bp, Bi, j, i, 0)] = BTx[j];
571                  }                  }
572                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);                  Free(rhs); Free(BTp); Free(BTx); Free(BTi);
# Line 730  Line 742
742
743      Memcpy(Bx, Ax, nnz);        /* x slot stays as is but w/o dim attribute */      Memcpy(Bx, Ax, nnz);        /* x slot stays as is but w/o dim attribute */
744
745      bdims[1] = ncb * adims[1];  /* bdims - need to find number of row blocks */      bdims[1] = ncb * adims[1];
746        /* find number of row blocks */
747      for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];      for (j = 0, nrb = -1; j < adims[2]; j++) if (Ai[j] > nrb) nrb = Ai[j];
748      bdims[0] = (nrb + 1) * adims[0]; /* +1 because of 0-based indices */      bdims[0] = (nrb + 1) * adims[0]; /* +1 because of 0-based indices */
749

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