[Distr-commits] r114 - pkg/distrMod/R

[noreply at r-forge.r-project.org]

Author: ruckdeschel
Date: 2008-04-04 13:31:42 +0200 (Fri, 04 Apr 2008)
New Revision: 114

Modified:
   pkg/distrMod/R/solve.R
Log:
minor enhancements to solve.R

Modified: pkg/distrMod/R/solve.R
===================================================================
--- pkg/distrMod/R/solve.R	2008-04-04 08:27:34 UTC (rev 113)
+++ pkg/distrMod/R/solve.R	2008-04-04 11:31:42 UTC (rev 114)
@@ -1,35 +1,27 @@
 setMethod("solve", signature(a = "ANY", b = "ANY"), function(a,b, generalized = TRUE,
-          tol = .Machine$double.eps, ...) {if(!generalized) return(base::solve(a,b, tol = tol, ...))
-          else if(is(try(return(base::solve(a,b, tol = tol, ...)), silent = TRUE),
-                                        "try-error")){
+          tol = .Machine$double.eps, ...) {
+                 if(!generalized) return(base::solve(a,b, tol = tol, ...))
+                 else if(is(try(return(base::solve(a,b, tol = tol, ...)), 
+                                silent = TRUE), "try-error")){
              if (!missing(b))
                 if(!(length(b)==nrow(a))) stop("non-conformable arguments")
-             a.svd <- svd(a)
-             d1 <- a.svd$d
-             d1.0 <- (d1 < tol) + 0.0
-             d1.1 <- 1/pmax(d1, d1.0)
-             d <- (1-d1.0) * d1.1
-             d <- if (length(d) == 1) d else diag(d)
-             a.m <- a.svd$v %*% d %*% t(a.svd$u)
-             erg <- if (missing(b)) a.m else a.m %*% b
-             return(erg)}
-             }
-             )
+             a.m <- ginv(a)
+             if (missing(b)) return(a.m) 
+             else return(a.m %*% b)
+             }})
 
-setMethod("solve", signature(a = "PosSemDefSymmMatrix", b = "ANY"), function(a,b, generalized = TRUE, tol = .Machine$double.eps, ...){
+setMethod("solve", signature(a = "PosSemDefSymmMatrix", b = "ANY"), 
+           function(a,b, generalized = TRUE, tol = .Machine$double.eps, ...){
           if(!generalized) return(base::solve(a,b, tol = tol, ...))
           else{
             er <- eigen(a)
             d1 <- er$values
-            d1.0 <- (d1 < tol) + 0.0
-            d1.1 <- 1/pmax(d1, d1.0)
-            d <- (1-d1.0) * d1.1
-            d <- if (length(d) == 1) d else diag(d)
-            A <- er$vectors %*% d %*% t(er$vectors)
+            d <- 1/d1[d1 > tol]
+            ev <- er$vectors[,d1 > tol]
+            A <- if (length(d)) ev %*% (t(ev)*d) else 0*a
             if(missing(b)) return(A)
-            else return(A%*%b)}
+            else return(A%*%b)}   
 })
 
 setMethod("solve", signature(a = "PosDefSymmMatrix", b = "ANY"), function(a,b, ...){
 base::solve(a,b, ...)})
-