Initial commit

Signed-off-by: Serguei Sokol <sokol@insa-toulouse.fr>

DESCRIPTION

 Package: gmresls
-Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9000
+Title: Solve Least Squares with GMRES(k)
+Version: 0.1.0
 Authors@R: 
-    person("First", "Last", , "first.last@example.com", role = c("aut", "cre"),
-           comment = c(ORCID = "YOUR-ORCID-ID"))
-Description: What the package does (one paragraph).
+    person("Serguei", "Sokol", , "sokol@insa-toulouse.fr", role = c("aut", "cre"),
+           comment = c(ORCID = "0000-0002-5674-3327"))
+Description: Solves a least squares system Ax=b (dim(A)=(m,n) with m >= n) with a preconditioner B: BAx=Bb (dim(B)=(n,m)). Implemented method uses GMRES(k) with callback functions, i.e. no explicit A or B are required. GMRES can be restarted after k iterations.
 License: GPL (>= 3)
 Encoding: UTF-8
 Roxygen: list(markdown = TRUE)
 RoxygenNote: 7.3.1
+Suggests: 
+    RUnit

NAMESPACE

 # Generated by roxygen2: do not edit by hand
 
+export(gmresls)

R/gmresls.R

+#' Solve a Least Squares System with a Preconditioner.
+#'
+#' Solve a least squares system Ax=b (dim(A)=(m,n) with m >= n) with a
+#' preconditioner B: BAx=Bb (dim(B)=(n,m)).
+#' Implemented method uses GMRES(k) with callback functions, i.e. no explicit A or B are required. GMRES can be restarted after k iterations.
+#'
+#' @param f_Ax A function calculating matrix-vector product Ax for a given x
+#' @param f_Br A function calculating matrix-vector product Br for a given r
+#' @param b A vector of right hand side
+#' @param x0 A vector or NULL (which means 0), initial approximation for Ax=b
+#' @param k An integer, parameter for restarting GMRES. Value 0 (default) means no restart, i.e. at most length(x) basis vectors will be constructed and used.
+#' @param maxit A maximal iteration number. Default (0) means length(x).
+#' @param tol A tolerance for solution, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7
+#' @param ... Parameters passed through to f_Ax and f_Br
+#' @returns The solution x, having the structure of Bb.
+#' @details
+#' Implemented method is equivalent to a classical GMRES(k) method with restart after constructing k basis vectors and applied to a square system BAx=Bb.
+#' Dense matrices constructed and stored by this method are of size (length(x), k) and (k+1, k) where k is GMRES current basis vector number. If maxit > k, GMRES will be restarted after each k iterations
+#' Particularity of this implementation that matrices A and B have no to be stored explicitly.
+#' User provides just callback function mimicking their multiplication by adequate vectors.
+#' In case of non convergence after maxit iterations, attr(x) will contain a field 'warning' with the message which will be also issued with warning()
+#' If the operator BA is not of full rank, iterations will be stopped before reaching convergence or maxit. A warning will be emitted in this case.
+#'
+#' @examples
+#' # prepare a 4x3 toy problem Ax=b
+#' A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
+#' xsol=1:3
+#' b=A%*%xsol
+#' f_Ax=function(x,...) list(...)$A%*%x
+#' f_Br=function(r,...) t(list(...)$A)%*%r
+#' x=gmresls(f_Ax, f_Br, b, A=A)
+#' @export
+gmresls=function(f_Ax, f_Br, b, x0=NULL, k=0, maxit=0, tol=1.e-7, ...) {
+  f_ba=function(x, ...) f_Br(f_Ax(x, ...), ...)
+  bb=f_Br(b, ...)
+  nbb=sqrt(sum(bb*bb)) # useful for convergence detection
+  nx=length(bb)
+  if (maxit == 0L)
+    maxit=nx
+  it=0L
+  warn=""
+  repeat { # restart till convergence or maxit reached
+    it=it+1L
+    r0=bb
+    if (!is.null(x0))
+      r0=r0-f_ba(x0, ...)
+    nr0=sqrt(sum(r0*r0))
+    v=as.matrix(r0/nr0)
+    h=matrix(0., 1L, 0L) # Heisenberg matrix, will be growing
+    re1=nr0
+    if (k == 0)
+      k=nx
+    converged=nr0/nbb <= tol
+    for (i in seq_len(k)) {
+      if (converged)
+        return(if (is.null(x0)) double(nx) else x0)
+      bavi=f_ba(v[,i], ...)
+      h=cbind(h, crossprod(v, bavi))
+      v=cbind(v, bavi-v%*%h[,i])
+      h=rbind(h, 0.)
+      h[i+1L,i]=sqrt(crossprod(v[,i+1]))
+      if (h[i+1L,i] <= tol) {
+        re1=c(re1, 0.)
+        y=qr.solve(h, re1)
+        nri=sqrt(crossprod(h%*%y-re1))
+        converged=nri/nr0 <= tol
+        if (it < nx)
+          warn="BA is not of full rank\n"
+        break
+      } else {
+        v[,i+1L]=v[,i+1L]/h[i+1L,i]
+      }
+      re1=c(re1, 0.)
+      y=qr.solve(h, re1)
+      nri=sqrt(crossprod(h%*%y-re1))
+      converged=nri/nr0 <= tol
+      it=it+1L
+    }
+    x=v[,-ncol(v)]%*%y
+    if( !is.null(x0) )
+      x=x+x0
+    # go to restart
+    if (converged || it >= maxit)
+      break
+  }
+  if (!converged)
+    warn=paste0(warn, "Maximal iteration number (", maxit, ") is reached without convergence")
+  if (nchar(warn))
+    warning(attr(x, "warning") <- warn)
+
+  attr(x, "resid")=nri
+  attr(x, "iteration")=it
+  x
+}

inst/unitTests/runit.gmresls.R

+# doc example
+# prepare a 4x3 toy problem Ax=b
+A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
+xsol=1:3
+b=A%*%xsol
+f_Ax=function(x,...) list(...)$A%*%x
+f_Br=function(r,...) t(list(...)$A)%*%r
+x=gmresls(f_Ax, f_Br, b, A=A)
+test.doc <- function() {
+  checkEqualsNumeric(x, xsol)
+}

man/gmresls.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/gmresls.R
+\name{gmresls}
+\alias{gmresls}
+\title{Solve a Least Squares System with a Preconditioner.}
+\usage{
+gmresls(f_Ax, f_Br, b, x0 = NULL, k = 0, maxit = 0, tol = 1e-07, ...)
+}
+\arguments{
+\item{f_Ax}{A function calculating matrix-vector product Ax for a given x}
+
+\item{f_Br}{A function calculating matrix-vector product Br for a given r}
+
+\item{b}{A vector of right hand side}
+
+\item{x0}{A vector or NULL (which means 0), initial approximation for Ax=b}
+
+\item{k}{An integer, parameter for restarting GMRES. Value 0 (default) means no restart, i.e. at most length(x) basis vectors will be constructed and used.}
+
+\item{maxit}{A maximal iteration number. Default (0) means length(x).}
+
+\item{tol}{A tolerance for solution, estimated as ||B(Ax-b)||/||Bb||, default 1.e-7}
+
+\item{...}{Parameters passed through to f_ax and f_br}
+}
+\value{
+The solution x, having the structure of Bb.
+}
+\description{
+Solve a least squares system Ax=b (dim(A)=(m,n) with m >= n) with a
+preconditioner B: BAx=Bb (dim(B)=(n,m)).
+Implemented method uses GMRES(k) with callback functions, i.e. no explicit A or B are required. GMRES can be restarted after k iterations.
+}
+\details{
+Implemented method is equivalent to a classical GMRES(k) method with restart after constructing k basis vectors and applied to a system BAx=Bb.
+Dense matrices constructed and stored by this method are of size (length(x), k) and (k+1, k) where k is GMRES current basis vector number. If maxit > k, GMRES will be restarted after each k iterations
+Particularity of this implementation that matrices A and B have no to be stored explicitly.
+User provides just callback function mimicking their multiplication by adequate vectors.
+In case of non convergence after maxit iterations, attr(x) will contain a field 'warning' with the message which will be also issued with warning()
+}
+\examples{
+# prepare a 4x3 toy problem Ax=b
+A=rbind(diag(1:3)+matrix(1, 3,3), rep(1, 3))
+xsol=1:3
+b=A\%*\%xsol
+f_Ax=function(x,...) list(...)$A\%*\%x
+f_Br=function(r,...) t(list(...)$A)\%*\%r
+x=gmresls(f_Ax, f_Br, b, A=A)
+}

tests/RUnit.R

+if (requireNamespace("RUnit", quietly=TRUE) && requireNamespace("gmresls", quietly=TRUE)) {
+   library(RUnit)
+   library(gmresls)
+
+   testSuite <- defineTestSuite(
+      name = "gmresls unit tests",
+      dirs = system.file("unitTests", package = "gmresls"),
+      testFuncRegexp = "^[Tt]est.+",
+      rngKind = RNGkind()[1L]
+   )
+   Sys.setenv("R_TESTS"="")
+   tests <- runTestSuite(testSuite)
+
+   printTextProtocol(tests)
+
+   if (getErrors(tests)$nFail > 0) stop("RUnit test failure")
+   if (getErrors(tests)$nErr > 0) stop("Errors in RUnit tests")
+}