diff --git a/DESCRIPTION b/DESCRIPTION
index 95df1379604ed2819fbe45a6bd6a333dd9718259..85bb15c8271d0b3885f20ca2a271575793ec2575 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,11 +1,13 @@
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
diff --git a/NAMESPACE b/NAMESPACE
index 6ae926839dd1829f1016a96f766d970ff184ad97..8396a214fd9faf09ab42ad04af3e4dc7fd651bcb 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,2 +1,3 @@
# Generated by roxygen2: do not edit by hand
+export(gmresls)
diff --git a/R/gmresls.R b/R/gmresls.R
new file mode 100644
index 0000000000000000000000000000000000000000..3f28766b6a8d039f14d46ab45002399159629064
--- /dev/null
+++ b/R/gmresls.R
@@ -0,0 +1,94 @@
+#' 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
+}
diff --git a/inst/unitTests/runit.gmresls.R b/inst/unitTests/runit.gmresls.R
new file mode 100644
index 0000000000000000000000000000000000000000..0289ff1ac0583127d01adc559f4cf7a68c44a031
--- /dev/null
+++ b/inst/unitTests/runit.gmresls.R
@@ -0,0 +1,11 @@
+# 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)
+}
diff --git a/man/gmresls.Rd b/man/gmresls.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..36301c6ccc29a9cdd1527fcc9ae896a5f1014bd9
--- /dev/null
+++ b/man/gmresls.Rd
@@ -0,0 +1,49 @@
+% 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)
+}
diff --git a/tests/RUnit.R b/tests/RUnit.R
new file mode 100644
index 0000000000000000000000000000000000000000..f4ddaa69752e2b1a304b69036c14628011ffe836
--- /dev/null
+++ b/tests/RUnit.R
@@ -0,0 +1,18 @@
+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")
+}