SCM

SCM Repository

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

Diff of /pkg/src/dgeMatrix.c

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

revision 10, Mon Mar 22 20:20:05 2004 UTC revision 301, Sun Oct 17 05:06:04 2004 UTC
# Line 256  Line 256 
256      UNPROTECT(1);      UNPROTECT(1);
257      return val;      return val;
258  }  }
259    
260    SEXP geMatrix_svd(SEXP x, SEXP nnu, SEXP nnv)
261    {
262        int nu = asInteger(nnu),
263            nv = asInteger(nnv),
264            *dims = INTEGER(GET_SLOT(x, Matrix_DimSym));
265        double *xx = REAL(GET_SLOT(x, Matrix_xSym));
266        SEXP val = PROTECT(allocVector(VECSXP, 3));
267    
268        if (dims[0] && dims[1]) {
269            int m = dims[0], n = dims[1], mm = (m < n)?m:n,
270                lwork = -1, info;
271            int *iwork = Calloc(8 * mm, int);
272            double tmp, *work;
273            int bdspac = 3*m*m + 4*m,
274                wrkbl, maxwrk, minwrk, itmp,
275                ione = 1, iminus1 = -1;
276            int i1, i2, i3;
277    
278            SET_VECTOR_ELT(val, 0, allocVector(REALSXP, mm));
279            SET_VECTOR_ELT(val, 1, allocMatrix(REALSXP, m, mm));
280            SET_VECTOR_ELT(val, 2, allocMatrix(REALSXP, mm, n));
281            F77_CALL(dgesdd)("S", &m, &n, xx, &m,
282                             REAL(VECTOR_ELT(val, 0)),
283                             REAL(VECTOR_ELT(val, 1)), &m,
284                             REAL(VECTOR_ELT(val, 2)), &mm,
285                             &tmp, &lwork, iwork, &info);
286            lwork = (int) tmp;
287            F77_CALL(foo)(&i1, &i2, &i3);
288            wrkbl = 3*m+(m+n)*i1;
289            if (wrkbl < (itmp = 3*m + m*i2)) wrkbl = itmp;
290            if (wrkbl < (itmp = 3*m + m*i3)) wrkbl = itmp;
291            itmp = bdspac+3*m;
292            maxwrk = (wrkbl > itmp) ? wrkbl : itmp;
293            minwrk = 3*m + ((bdspac > n) ?  bdspac : n);
294            work = Calloc(lwork, double);
295            F77_CALL(dgesdd)("S", &m, &n, xx, &m,
296                             REAL(VECTOR_ELT(val, 0)),
297                             REAL(VECTOR_ELT(val, 1)), &m,
298                             REAL(VECTOR_ELT(val, 2)), &mm,
299                             work, &lwork, iwork, &info);
300            Free(iwork); Free(work);
301        }
302        UNPROTECT(1);
303        return val;
304    }

Legend:
Removed from v.10  
changed lines
  Added in v.301

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