From 20a2ca709f107e9e5e8c66fde554907f27517799 Mon Sep 17 00:00:00 2001 From: jacobwilliams Date: Thu, 25 Jan 2024 04:21:44 +0000 Subject: [PATCH] =?UTF-8?q?Deploying=20to=20gh-pages=20from=20@=20jacobwil?= =?UTF-8?q?liams/LSMR@7364adbcc4755580c49d448d2b41875cead7745b=20?= =?UTF-8?q?=F0=9F=9A=80?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- module/lsmrblas.html | 5 + modules.json | 2 +- sourcefile/lsmrblas.f90.html | 487 +++++++++++++++-------------- src/lsmrblas.F90 | 9 +- tipuesearch/tipuesearch_content.js | 2 +- 5 files changed, 262 insertions(+), 243 deletions(-) diff --git a/module/lsmrblas.html b/module/lsmrblas.html index b1e6442..a0f498c 100644 --- a/module/lsmrblas.html +++ b/module/lsmrblas.html @@ -172,6 +172,11 @@

References

ACM Transactions on Mathematical Software, Volume 5, Number 3, September 1979, pages 308-323. +
+

Note

+

This version has not been linked with an external BLAS library, + and uses the BLAS functions in this module.

+

Uses

diff --git a/modules.json b/modules.json index 2da4568..a101500 100644 --- a/modules.json +++ b/modules.json @@ -1 +1 @@ -{"ford-metadata": {"version": "7.0.5"}, "modules": [{"name": "lsmrDataModule", "external_url": "./module/lsmrdatamodule.html", "obj": "module", "pub_procs": {}, "pub_absints": {}, "pub_types": {}, "pub_vars": {"lsmr_wp": {"name": "lsmr_wp", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_wp", "obj": "variable", "vartype": "integer", "permission": "public"}, "lsmr_ip": {"name": "lsmr_ip", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_ip", "obj": "variable", "vartype": "integer", "permission": "public"}}, "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [{"name": "lsmr_wp", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_wp", "obj": "variable", "vartype": "integer", "permission": "public"}, {"name": "lsmr_ip", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_ip", "obj": "variable", "vartype": "integer", "permission": "public"}], "permission": "private"}, {"name": "lsmrModule", "external_url": "./module/lsmrmodule.html", "obj": "module", "pub_procs": {"lsmr": {"name": "lsmr", "external_url": "./proc/lsmr.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "lsmr_ez": {"name": "lsmr_ez", "external_url": "./proc/lsmr_ez.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}}, "pub_absints": {}, "pub_types": {}, "pub_vars": {}, "functions": [], "subroutines": [{"name": "lsmr", "external_url": "./proc/lsmr.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "lsmr_ez", "external_url": "./proc/lsmr_ez.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "interfaces": [], "absinterfaces": [{"name": "Aprod1_f", "external_url": "./interface/aprod1_f.html", "obj": "interface", "proctype": "Interface", "variables": [], "permission": "private", "generic": "False"}, {"name": "Aprod2_f", "external_url": "./interface/aprod2_f.html", "obj": "interface", "proctype": "Interface", "variables": [], "permission": "private", "generic": "False"}], "types": [], "variables": [], "permission": "private"}, {"name": "lsmrblas", "external_url": "./module/lsmrblas.html", "obj": "module", "pub_procs": {"ddot": {"name": "ddot", "external_url": "./proc/ddot.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dnrm2": {"name": "dnrm2", "external_url": "./proc/dnrm2.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dcopy": {"name": "dcopy", "external_url": "./proc/dcopy.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dscal": {"name": "dscal", "external_url": "./proc/dscal.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}}, "pub_absints": {}, "pub_types": {}, "pub_vars": {}, "functions": [{"name": "ddot", "external_url": "./proc/ddot.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "dnrm2", "external_url": "./proc/dnrm2.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "subroutines": [{"name": "dcopy", "external_url": "./proc/dcopy.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "dscal", "external_url": "./proc/dscal.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "private"}]} \ No newline at end of file +{"ford-metadata": {"version": "7.0.5"}, "modules": [{"name": "lsmrModule", "external_url": "./module/lsmrmodule.html", "obj": "module", "pub_procs": {"lsmr": {"name": "lsmr", "external_url": "./proc/lsmr.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "lsmr_ez": {"name": "lsmr_ez", "external_url": "./proc/lsmr_ez.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}}, "pub_absints": {}, "pub_types": {}, "pub_vars": {}, "functions": [], "subroutines": [{"name": "lsmr", "external_url": "./proc/lsmr.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "lsmr_ez", "external_url": "./proc/lsmr_ez.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "interfaces": [], "absinterfaces": [{"name": "Aprod1_f", "external_url": "./interface/aprod1_f.html", "obj": "interface", "proctype": "Interface", "variables": [], "permission": "private", "generic": "False"}, {"name": "Aprod2_f", "external_url": "./interface/aprod2_f.html", "obj": "interface", "proctype": "Interface", "variables": [], "permission": "private", "generic": "False"}], "types": [], "variables": [], "permission": "private"}, {"name": "lsmrblas", "external_url": "./module/lsmrblas.html", "obj": "module", "pub_procs": {"ddot": {"name": "ddot", "external_url": "./proc/ddot.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dnrm2": {"name": "dnrm2", "external_url": "./proc/dnrm2.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dcopy": {"name": "dcopy", "external_url": "./proc/dcopy.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, "dscal": {"name": "dscal", "external_url": "./proc/dscal.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}}, "pub_absints": {}, "pub_types": {}, "pub_vars": {}, "functions": [{"name": "ddot", "external_url": "./proc/ddot.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "dnrm2", "external_url": "./proc/dnrm2.html", "obj": "proc", "proctype": "Function", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "subroutines": [{"name": "dcopy", "external_url": "./proc/dcopy.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}, {"name": "dscal", "external_url": "./proc/dscal.html", "obj": "proc", "proctype": "Subroutine", "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "public"}], "interfaces": [], "absinterfaces": [], "types": [], "variables": [], "permission": "private"}, {"name": "lsmrDataModule", "external_url": "./module/lsmrdatamodule.html", "obj": "module", "pub_procs": {}, "pub_absints": {}, "pub_types": {}, "pub_vars": {"lsmr_wp": {"name": "lsmr_wp", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_wp", "obj": "variable", "vartype": "integer", "permission": "public"}, "lsmr_ip": {"name": "lsmr_ip", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_ip", "obj": "variable", "vartype": "integer", "permission": "public"}}, "functions": [], "subroutines": [], "interfaces": [], "absinterfaces": [], "types": [], "variables": [{"name": "lsmr_wp", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_wp", "obj": "variable", "vartype": "integer", "permission": "public"}, {"name": "lsmr_ip", "external_url": "./module/lsmrdatamodule.html#variable-lsmr_ip", "obj": "variable", "vartype": "integer", "permission": "public"}], "permission": "private"}]} \ No newline at end of file diff --git a/sourcefile/lsmrblas.f90.html b/sourcefile/lsmrblas.f90.html index 4d8a9e7..58b8425 100644 --- a/sourcefile/lsmrblas.f90.html +++ b/sourcefile/lsmrblas.f90.html @@ -256,258 +256,265 @@

Source Code

! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. - - module lsmrblas - - use lsmrDataModule, only : ip => lsmr_ip, wp => lsmr_wp - - implicit none - - private +! +#ifdef HAS_BLAS +!@note This version has been linked with an external BLAS library, +! only the interfaces are defined in this module. +#else +!@note This version has not been linked with an external BLAS library, +! and uses the BLAS functions in this module. +#endif - public :: ddot, dnrm2, dscal, dcopy + module lsmrblas -#ifdef HAS_BLAS + use lsmrDataModule, only : ip => lsmr_ip, wp => lsmr_wp - ! get BLAS from an externally-linked library (double precision only) + implicit none - interface + private - function ddot (n,dx,incx,dy,incy) - implicit none - integer, intent(in) :: n,incx,incy - double precision, intent(in) :: dx(*),dy(*) - double precision :: ddot - end function ddot - - function dnrm2 (n,dx,incx) - implicit none - integer, intent(in) :: n,incx - double precision, intent(in) :: dx(*) - double precision :: dnrm2 - end function dnrm2 - - subroutine dscal (n,sa,x,incx) - implicit none - integer, intent(in) :: n,incx - double precision, intent(in) :: sa - double precision, intent(inout) :: x(*) - end subroutine dscal - - subroutine dcopy(n,dx,incx,dy,incy) - implicit none - double precision dx(*),dy(*) - integer i,incx,incy,ix,iy,m,n - end subroutine dcopy - - end interface + public :: ddot, dnrm2, dscal, dcopy + +#ifdef HAS_BLAS + + ! get BLAS from an externally-linked library (double precision only) + + interface + + function ddot (n,dx,incx,dy,incy) + implicit none + integer, intent(in) :: n,incx,incy + double precision, intent(in) :: dx(*),dy(*) + double precision :: ddot + end function ddot + + function dnrm2 (n,dx,incx) + implicit none + integer, intent(in) :: n,incx + double precision, intent(in) :: dx(*) + double precision :: dnrm2 + end function dnrm2 + + subroutine dscal (n,sa,x,incx) + implicit none + integer, intent(in) :: n,incx + double precision, intent(in) :: sa + double precision, intent(inout) :: x(*) + end subroutine dscal - -#else - - contains - -!***************************************************************************** -!> -! Copies a vector X to a vector Y. - - subroutine dcopy(n,dx,incx,dy,incy) - - implicit none - real(wp) dx(*),dy(*) - integer(ip) i,incx,incy,ix,iy,m,n - - if ( n <= 0 ) then - return - end if - - if ( incx == 1 .and. incy == 1 ) then - - m = mod ( n, 7 ) - - if ( m /= 0 ) then - dy(1:m) = dx(1:m) - end if - - do i = m+1, n, 7 - dy(i) = dx(i) - dy(i + 1) = dx(i + 1) - dy(i + 2) = dx(i + 2) - dy(i + 3) = dx(i + 3) - dy(i + 4) = dx(i + 4) - dy(i + 5) = dx(i + 5) - dy(i + 6) = dx(i + 6) - end do - - else - - if ( 0 <= incx ) then - ix = 1 - else - ix = ( -n + 1 ) * incx + 1 - end if - - if ( 0 <= incy ) then - iy = 1 - else - iy = ( -n + 1 ) * incy + 1 - end if - - do i = 1, n - dy(iy) = dx(ix) - ix = ix + incx - iy = iy + incy - end do - end if + subroutine dcopy(n,dx,incx,dy,incy) + implicit none + double precision dx(*),dy(*) + integer i,incx,incy,ix,iy,m,n + end subroutine dcopy + + end interface + +#else + + contains + +!***************************************************************************** +!> +! Copies a vector X to a vector Y. + + subroutine dcopy(n,dx,incx,dy,incy) + + implicit none + real(wp) dx(*),dy(*) + integer(ip) i,incx,incy,ix,iy,m,n + + if ( n <= 0 ) then + return + end if + + if ( incx == 1 .and. incy == 1 ) then + + m = mod ( n, 7 ) + + if ( m /= 0 ) then + dy(1:m) = dx(1:m) + end if + + do i = m+1, n, 7 + dy(i) = dx(i) + dy(i + 1) = dx(i + 1) + dy(i + 2) = dx(i + 2) + dy(i + 3) = dx(i + 3) + dy(i + 4) = dx(i + 4) + dy(i + 5) = dx(i + 5) + dy(i + 6) = dx(i + 6) + end do + + else + + if ( 0 <= incx ) then + ix = 1 + else + ix = ( -n + 1 ) * incx + 1 + end if + + if ( 0 <= incy ) then + iy = 1 + else + iy = ( -n + 1 ) * incy + 1 + end if -end subroutine dcopy - -!***************************************************************************** -!> -! Dot product of two vectors. - - real(wp) function ddot(n,dx,incx,dy,incy) - - integer(ip),intent(in) :: n !! the number of entries in the vectors. - real(wp),intent(in) :: dx(*) !! the first vector - integer(ip),intent(in) :: incx !! the increment between successive entries in DX. - real(wp),intent(in) :: dy(*) !! the second vector - integer(ip),intent(in) :: incy !! the increment between successive entries in DY. - - integer(ip) :: i,ix,iy,m - real(wp) :: dtemp - - ddot = 0.0_wp - dtemp = 0.0_wp - if ( n <= 0 ) return + do i = 1, n + dy(iy) = dx(ix) + ix = ix + incx + iy = iy + incy + end do + end if + +end subroutine dcopy + +!***************************************************************************** +!> +! Dot product of two vectors. + + real(wp) function ddot(n,dx,incx,dy,incy) + + integer(ip),intent(in) :: n !! the number of entries in the vectors. + real(wp),intent(in) :: dx(*) !! the first vector + integer(ip),intent(in) :: incx !! the increment between successive entries in DX. + real(wp),intent(in) :: dy(*) !! the second vector + integer(ip),intent(in) :: incy !! the increment between successive entries in DY. - if ( incx /= 1 .or. incy /= 1 ) then - ! Code for unequal increments or equal increments - ! not equal to 1. - - if ( 0 <= incx ) then - ix = 1 - else - ix = ( - n + 1 ) * incx + 1 - end if - - if ( 0 <= incy ) then - iy = 1 - else - iy = ( - n + 1 ) * incy + 1 - end if - - do i = 1, n - dtemp = dtemp + dx(ix) * dy(iy) - ix = ix + incx - iy = iy + incy - end do - - else - ! Code for both increments equal to 1. - - m = mod ( n, 5 ) - - do i = 1, m - dtemp = dtemp + dx(i) * dy(i) - end do - - do i = m+1, n, 5 - dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) & - + dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4) - end do - - end if + integer(ip) :: i,ix,iy,m + real(wp) :: dtemp + + ddot = 0.0_wp + dtemp = 0.0_wp + if ( n <= 0 ) return + + if ( incx /= 1 .or. incy /= 1 ) then + ! Code for unequal increments or equal increments + ! not equal to 1. + + if ( 0 <= incx ) then + ix = 1 + else + ix = ( - n + 1 ) * incx + 1 + end if + + if ( 0 <= incy ) then + iy = 1 + else + iy = ( - n + 1 ) * incy + 1 + end if + + do i = 1, n + dtemp = dtemp + dx(ix) * dy(iy) + ix = ix + incx + iy = iy + incy + end do + + else + ! Code for both increments equal to 1. + + m = mod ( n, 5 ) + + do i = 1, m + dtemp = dtemp + dx(i) * dy(i) + end do - ddot = dtemp - -end function ddot - -!***************************************************************************** -!> -! The euclidean norm of a vector `sqrt ( X' * X )`. - - real(wp) function dnrm2 ( n, x, incx) - - integer(ip),intent(in) :: n - real(wp),intent(in) :: x(*) - integer(ip),intent(in) :: incx - - integer(ip) :: ix - real(wp) :: ssq,absxi,norm,scale + do i = m+1, n, 5 + dtemp = dtemp + dx(i)*dy(i) + dx(i+1)*dy(i+1) + dx(i+2)*dy(i+2) & + + dx(i+3)*dy(i+3) + dx(i+4)*dy(i+4) + end do + + end if + + ddot = dtemp + +end function ddot + +!***************************************************************************** +!> +! The euclidean norm of a vector `sqrt ( X' * X )`. + + real(wp) function dnrm2 ( n, x, incx) - if ( n < 1 .or. incx < 1 ) then - norm = 0.0_wp - else if ( n == 1 ) then - norm = abs ( x(1) ) - else - scale = 0.0_wp - ssq = 1.0_wp - - do ix = 1, 1 + ( n - 1 )*incx, incx - if ( x(ix) /= 0.0_wp ) then - absxi = abs ( x(ix) ) - if ( scale < absxi ) then - ssq = 1.0_wp + ssq * ( scale / absxi )**2 - scale = absxi - else - ssq = ssq + ( absxi / scale )**2 - end if - end if - end do - norm = scale * sqrt ( ssq ) - end if - - dnrm2 = norm - -end function dnrm2 - -!***************************************************************************** -!> -! Scales a vector by a constant. - - subroutine dscal(n,sa,x,incx) - - integer(ip),intent(in) :: n !! the number of entries in the vector. - real(wp) ,intent(in) :: sa !! the multiplier. - real(wp),intent(inout) :: x(*) !! the vector to be scaled. - integer(ip),intent(in) :: incx !! the increment between successive entries of X. + integer(ip),intent(in) :: n + real(wp),intent(in) :: x(*) + integer(ip),intent(in) :: incx + + integer(ip) :: ix + real(wp) :: ssq,absxi,norm,scale + + if ( n < 1 .or. incx < 1 ) then + norm = 0.0_wp + else if ( n == 1 ) then + norm = abs ( x(1) ) + else + scale = 0.0_wp + ssq = 1.0_wp + + do ix = 1, 1 + ( n - 1 )*incx, incx + if ( x(ix) /= 0.0_wp ) then + absxi = abs ( x(ix) ) + if ( scale < absxi ) then + ssq = 1.0_wp + ssq * ( scale / absxi )**2 + scale = absxi + else + ssq = ssq + ( absxi / scale )**2 + end if + end if + end do + norm = scale * sqrt ( ssq ) + end if + + dnrm2 = norm + +end function dnrm2 + +!***************************************************************************** +!> +! Scales a vector by a constant. - integer(ip) :: i, ix - integer(ip) :: m - - if ( n <= 0 ) then - return - else if ( incx == 1 ) then - m = mod ( n, 5 ) - x(1:m) = sa * x(1:m) - - do i = m+1, n, 5 - x(i) = sa * x(i) - x(i+1) = sa * x(i+1) - x(i+2) = sa * x(i+2) - x(i+3) = sa * x(i+3) - x(i+4) = sa * x(i+4) - end do - else - if ( 0 <= incx ) then - ix = 1 - else - ix = ( - n + 1 ) * incx + 1 - end if - - do i = 1, n - x(ix) = sa * x(ix) - ix = ix + incx - end do - - end if + subroutine dscal(n,sa,x,incx) + + integer(ip),intent(in) :: n !! the number of entries in the vector. + real(wp) ,intent(in) :: sa !! the multiplier. + real(wp),intent(inout) :: x(*) !! the vector to be scaled. + integer(ip),intent(in) :: incx !! the increment between successive entries of X. + + integer(ip) :: i, ix + integer(ip) :: m + + if ( n <= 0 ) then + return + else if ( incx == 1 ) then + m = mod ( n, 5 ) + x(1:m) = sa * x(1:m) + + do i = m+1, n, 5 + x(i) = sa * x(i) + x(i+1) = sa * x(i+1) + x(i+2) = sa * x(i+2) + x(i+3) = sa * x(i+3) + x(i+4) = sa * x(i+4) + end do + else + if ( 0 <= incx ) then + ix = 1 + else + ix = ( - n + 1 ) * incx + 1 + end if -end subroutine dscal - -#endif - -end module lsmrblas + do i = 1, n + x(ix) = sa * x(ix) + ix = ix + incx + end do + + end if + +end subroutine dscal + +#endif + +end module lsmrblas
diff --git a/src/lsmrblas.F90 b/src/lsmrblas.F90 index a25fe40..9a25d39 100644 --- a/src/lsmrblas.F90 +++ b/src/lsmrblas.F90 @@ -16,6 +16,14 @@ ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. +! +#ifdef HAS_BLAS +!@note This version has been linked with an external BLAS library, +! only the interfaces are defined in this module. +#else +!@note This version has not been linked with an external BLAS library, +! and uses the BLAS functions in this module. +#endif module lsmrblas @@ -62,7 +70,6 @@ end subroutine dcopy end interface - #else contains diff --git a/tipuesearch/tipuesearch_content.js b/tipuesearch/tipuesearch_content.js index 32abd29..79bf834 100644 --- a/tipuesearch/tipuesearch_content.js +++ b/tipuesearch/tipuesearch_content.js @@ -1 +1 @@ -var tipuesearch = {"pages":[{"title":" LSMR ","text":"LSMR LSMR: Sparse Equations and Least Squares. This is a slightly modernized version of the original FSMR (f90 version) by David Fong & Michael Saunders, Systems Optimization Laboratory (SOL), Stanford University, Stanford, CA 94305-4026, USA. Status Compiling A Fortran Package Manager manifest file is included, so that the library and test cases can be compiled with FPM. For example: fpm build --profile release\nfpm test --profile release To use lsmr within your fpm project, add the following to your fpm.toml file: [dependencies] LSMR = { git = \"https://github.com/jacobwilliams/LSMR.git\" } License The original version of LSMR (f90 version) was provided by SOL, Stanford University under the terms of the OSI Common Public License (CPL) or the BSD License . Documentation The latest API documentation can be found here . This was generated from the source code using FORD . See also LSMR: Sparse Equations and Least Squares -- LSMR original code. Sparse Parallel Robust Algorithms Library -- another version Developer Info Jacob Williams","tags":"home","loc":"index.html"},{"title":"Aprod1_f – LSMR","text":"interface private subroutine Aprod1_f(m, n, x, y) Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(in) :: x (n) real(kind=dp), intent(inout) :: y (m) Description y := y + A*x","tags":"","loc":"interface/aprod1_f.html"},{"title":"Aprod2_f – LSMR","text":"interface private subroutine Aprod2_f(m, n, x, y) Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(inout) :: x (n) real(kind=dp), intent(in) :: y (m) Description x := x + A'*y","tags":"","loc":"interface/aprod2_f.html"},{"title":"lsmr – LSMR","text":"public subroutine lsmr(m, n, Aprod1, Aprod2, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) LSMR finds a solution x to the following problems: Unsymmetric equations: Solve A*x = b Linear least squares: Solve A*x = b in the least-squares sense Damped least squares: Solve ( A )*x = ( b )\n ( damp*I ) ( 0 ) in the least-squares sense where A is a matrix with m rows and n columns, b is an m-vector,\n and damp is a scalar. (All quantities are real.)\n The matrix A is treated as a linear operator. It is accessed\n by means of subroutine calls with the following purpose: call Aprod1(m,n,x,y) must compute y = y + A*x without altering x . call Aprod2(m,n,x,y) must compute x = x + A'*y without altering y . LSMR uses an iterative method to approximate the solution.\n The number of iterations required to reach a certain accuracy\n depends strongly on the scaling of the problem. Poor scaling of\n the rows or columns of A should therefore be avoided where\n possible. For example, in problem 1 the solution is unaltered by\n row-scaling. If a row of A is very small or large compared to\n the other rows of A, the corresponding row of ( A b ) should be\n scaled up or down. In problems 1 and 2, the solution x is easily recovered\n following column-scaling. Unless better information is known,\n the nonzero columns of A should be scaled so that they all have\n the same Euclidean norm (e.g., 1.0). In problem 3, there is no freedom to re-scale if damp is\n nonzero. However, the value of damp should be assigned only\n after attention has been paid to the scaling of A. The parameter damp is intended to help regularize\n ill-conditioned systems, by preventing the true solution from\n being very large. Another aid to regularization is provided by\n the parameter condA, which may be used to terminate iterations\n before the computed solution becomes very large. Note that x is not an input parameter.\n If some initial estimate x0 is known and if damp = 0,\n one could proceed as follows: Compute a residual vector r0 = b - A*x0 . Use LSMR to solve the system A*dx = r0 . Add the correction dx to obtain a final solution x = x0 + dx . This requires that x0 be available before and after the call\n to LSMR. To judge the benefits, suppose LSMR takes k1 iterations\n to solve A*x = b and k2 iterations to solve A*dx = r0 .\n If x0 is \"good\", norm(r0) will be smaller than norm(b).\n If the same stopping tolerances atol and btol are used for each\n system, k1 and k2 will be similar, but the final solution x0 + dx should be more accurate. The only way to reduce the total work\n is to use a larger stopping tolerance for the second system.\n If some value btol is suitable for A*x = b , the larger value btol*norm(b)/norm(r0) should be suitable for A*dx = r0 . Preconditioning is another way to reduce the number of iterations.\n If it is possible to solve a related system M*x = b efficiently,\n where M approximates A in some helpful way\n (e.g. M - A has low rank or its elements are small relative to\n those of A), LSMR may converge more rapidly on the system A*M(inverse)*z = b ,\n after which x can be recovered by solving M*x = z . NOTE: If A is symmetric, LSMR should not be used!\n Alternatives are the symmetric conjugate-gradient method (CG)\n and/or SYMMLQ.\n SYMMLQ is an implementation of symmetric CG that applies to\n any symmetric A and will converge more rapidly than LSMR.\n If A is positive definite, there are other implementations of\n symmetric CG that require slightly less work per iteration\n than SYMMLQ (but will take the same number of iterations). Notation The following quantities are used in discussing the subroutine\n parameters: Abar = ( A ), bbar = (b)\n (damp*I) (0)\n\n r = b - A*x, rbar = bbar - Abar*x\n\n normr = sqrt( norm(r)**2 + damp**2 * norm(x)* *2 )\n = norm( rbar )\n\n eps = the relative precision of floating-point arithmetic.\n On most machines, eps is about 1.0e-7 and 1.0e-16\n in single and double precision respectively.\n We expect eps to be about 1e-16 always. LSMR minimizes the function normr with respect to x . Precision The number of iterations required by LSMR will decrease\n if the computation is performed in higher precision. Reference http://www.stanford.edu/group/SOL/software/lsmr.html LSMR development: 21 Sep 2007: Fortran 90 version of LSQR implemented.\n Aprod1, Aprod2 implemented via f90 interface. 17 Jul 2010: LSMR derived from LSQR and lsmr.m. 07 Sep 2010: Local reorthogonalization now working. 02 May 2014: With damp>0, istop=2 was incorrectly set to istop=3\n (so incorrect stopping message was printed). Fixed. Note Any or all of atol , btol , conlim may be set to zero.\n The effect will be the same as the values eps , eps , 1/eps . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m the number of rows in A. integer(kind=ip), intent(in) :: n the number of columns in A. procedure( Aprod1_f ) :: Aprod1 See above. procedure( Aprod2_f ) :: Aprod2 See above. real(kind=dp), intent(in) :: b (m) The rhs vector b . real(kind=dp), intent(in) :: damp The damping parameter for problem 3 above.\n(damp should be 0.0 for problems 1 and 2.)\nIf the system A*x = b is incompatible, values\nof damp in the range 0 to sqrt(eps)*norm(A) will probably have a negligible effect.\nLarger values of damp will tend to decrease\nthe norm of x and reduce the number of\niterations required by LSMR. The work per iteration and the storage needed\nby LSMR are the same for all values of damp. real(kind=dp), intent(in) :: atol An estimate of the relative error in the data\ndefining the matrix A. For example, if A is\naccurate to about 6 digits, set atol = 1.0e-6. real(kind=dp), intent(in) :: btol An estimate of the relative error in the data\ndefining the rhs b. For example, if b is\naccurate to about 6 digits, set btol = 1.0e-6. real(kind=dp), intent(in) :: conlim An upper limit on cond(Abar), the apparent\ncondition number of the matrix Abar.\nIterations will be terminated if a computed\nestimate of cond(Abar) exceeds conlim.\nThis is intended to prevent certain small or\nzero singular values of A or Abar from\ncoming into effect and causing unwanted growth\nin the computed solution. conlim and damp may be used separately or\ntogether to regularize ill-conditioned systems. Normally, conlim should be in the range\n1000 to 1/eps.\nSuggested value: conlim = 1/(100*eps) for compatible systems, conlim = 1/(10*sqrt(eps)) for least squares. integer(kind=ip), intent(in) :: itnlim An upper limit on the number of iterations.\nSuggested value: itnlim = n/2 for well-conditioned systems\n with clustered singular values, itnlim = 4*n otherwise. integer(kind=ip), intent(in) :: localSize No. of vectors for local reorthogonalization: 0 No reorthogonalization is performed. 0 This many n-vectors \"v\" (the most recent ones)\n are saved for reorthogonalizing the next v. localSize need not be more than min(m,n).\nAt most min(m,n) vectors will be allocated. integer(kind=ip), intent(in) :: nout File number for printed output. If positive,\na summary will be printed on file nout. real(kind=dp), intent(out) :: x (n) Returns the computed solution x . integer(kind=ip), intent(out) :: istop An integer giving the reason for termination: 0 x = 0 is the exact solution.\n No iterations were performed. 1 The equations A*x = b are probably compatible. Norm(A*x - b) is sufficiently small, given the\n values of atol and btol. 2 damp is zero. The system A*x = b is probably\n not compatible. A least-squares solution has\n been obtained that is sufficiently accurate,\n given the value of atol. 3 damp is nonzero. A damped least-squares\n solution has been obtained that is sufficiently\n accurate, given the value of atol. 4 An estimate of cond(Abar) has exceeded conlim.\n The system A*x = b appears to be ill-conditioned,\n or there could be an error in Aprod1 or Aprod2. 5 The iteration limit itnlim was reached. integer(kind=ip), intent(out) :: itn The number of iterations performed. real(kind=dp), intent(out) :: normA An estimate of the Frobenius norm of Abar.\nThis is the square-root of the sum of squares\nof the elements of Abar.\nIf damp is small and the columns of A\nhave all been scaled to have length 1.0,\nnormA should increase to roughly sqrt(n).\nA radically different value for normA may\nindicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: condA An estimate of cond(Abar), the condition\nnumber of Abar. A very high value of condA\nmay again indicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: normr An estimate of the final value of norm(rbar),\nthe function being minimized (see notation\nabove). This will be small if A*x = b has\na solution. real(kind=dp), intent(out) :: normAr An estimate of the final value of norm( Abar'*rbar ) , the norm of\nthe residual for the normal equations.\nThis should be small in all cases. (normAr\nwill often be smaller than the true value\ncomputed from the output vector x.) real(kind=dp), intent(out) :: normx An estimate of norm(x) for the final solution x. Called by proc~~lsmr~~CalledByGraph proc~lsmr lsmrModule::lsmr proc~lsmr_ez lsmrModule::lsmr_ez proc~lsmr_ez->proc~lsmr Help Graph Key Nodes of different colours represent the following: Graph Key Subroutine Subroutine Function Function Interface Interface Type Bound Procedure Type Bound Procedure Unknown Procedure Type Unknown Procedure Type Program Program This Page's Entity This Page's Entity Solid arrows point from a procedure to one which it calls. Dashed \narrows point from an interface to procedures which implement that interface.\nThis could include the module procedures in a generic interface or the\nimplementation in a submodule of an interface in a parent module. Source Code subroutine lsmr ( m , n , Aprod1 , Aprod2 , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer ( ip ), intent ( in ) :: m !! the number of rows in A. integer ( ip ), intent ( in ) :: n !! the number of columns in A. integer ( ip ), intent ( in ) :: itnlim !! An upper limit on the number of iterations. !! Suggested value: !! !! * `itnlim = n/2` for well-conditioned systems !! with clustered singular values, !! * `itnlim = 4*n` otherwise. integer ( ip ), intent ( in ) :: localSize !! No. of vectors for local reorthogonalization: !! !! * 0 No reorthogonalization is performed. !! * >0 This many n-vectors \"v\" (the most recent ones) !! are saved for reorthogonalizing the next v. !! !! localSize need not be more than min(m,n). !! At most min(m,n) vectors will be allocated. integer ( ip ), intent ( in ) :: nout !! File number for printed output. If positive, !! a summary will be printed on file nout. integer ( ip ), intent ( out ) :: istop !! An integer giving the reason for termination: !! !! * 0 x = 0 is the exact solution. !! No iterations were performed. !! * 1 The equations `A*x = b` are probably compatible. !! `Norm(A*x - b)` is sufficiently small, given the !! values of atol and btol. !! * 2 damp is zero. The system `A*x = b` is probably !! not compatible. A least-squares solution has !! been obtained that is sufficiently accurate, !! given the value of atol. !! * 3 damp is nonzero. A damped least-squares !! solution has been obtained that is sufficiently !! accurate, given the value of atol. !! * 4 An estimate of cond(Abar) has exceeded conlim. !! The system `A*x = b` appears to be ill-conditioned, !! or there could be an error in Aprod1 or Aprod2. !! * 5 The iteration limit itnlim was reached. integer ( ip ), intent ( out ) :: itn !! The number of iterations performed. real ( dp ), intent ( in ) :: b ( m ) !! The rhs vector `b`. real ( dp ), intent ( out ) :: x ( n ) !! Returns the computed solution `x`. real ( dp ), intent ( in ) :: atol !! An estimate of the relative error in the data !! defining the matrix A. For example, if A is !! accurate to about 6 digits, set atol = 1.0e-6. real ( dp ), intent ( in ) :: btol !! An estimate of the relative error in the data !! defining the rhs b. For example, if b is !! accurate to about 6 digits, set btol = 1.0e-6. real ( dp ), intent ( in ) :: conlim !! An upper limit on cond(Abar), the apparent !! condition number of the matrix Abar. !! Iterations will be terminated if a computed !! estimate of cond(Abar) exceeds conlim. !! This is intended to prevent certain small or !! zero singular values of A or Abar from !! coming into effect and causing unwanted growth !! in the computed solution. !! !! conlim and damp may be used separately or !! together to regularize ill-conditioned systems. !! !! Normally, conlim should be in the range !! 1000 to 1/eps. !! Suggested value: !! !! * `conlim = 1/(100*eps)` for compatible systems, !! * `conlim = 1/(10*sqrt(eps))` for least squares. real ( dp ), intent ( in ) :: damp !! The damping parameter for problem 3 above. !! (damp should be 0.0 for problems 1 and 2.) !! If the system `A*x = b` is incompatible, values !! of damp in the range 0 to `sqrt(eps)*norm(A)` !! will probably have a negligible effect. !! Larger values of damp will tend to decrease !! the norm of x and reduce the number of !! iterations required by LSMR. !! !! The work per iteration and the storage needed !! by LSMR are the same for all values of damp. real ( dp ), intent ( out ) :: normA !! An estimate of the Frobenius norm of Abar. !! This is the square-root of the sum of squares !! of the elements of Abar. !! If damp is small and the columns of A !! have all been scaled to have length 1.0, !! normA should increase to roughly sqrt(n). !! A radically different value for normA may !! indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: condA !! An estimate of cond(Abar), the condition !! number of Abar. A very high value of condA !! may again indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: normr !! An estimate of the final value of norm(rbar), !! the function being minimized (see notation !! above). This will be small if A*x = b has !! a solution. real ( dp ), intent ( out ) :: normAr !! An estimate of the final value of !! `norm( Abar'*rbar )`, the norm of !! the residual for the normal equations. !! This should be small in all cases. (normAr !! will often be smaller than the true value !! computed from the output vector x.) real ( dp ), intent ( out ) :: normx !! An estimate of norm(x) for the final solution x. procedure ( Aprod1_f ) :: Aprod1 !! See above. procedure ( Aprod2_f ) :: Aprod2 !! See above. ! Local arrays and variables real ( dp ) :: h ( n ), hbar ( n ), u ( m ), v ( n ), w ( n ), localV ( n , min ( localSize , m , n )) logical :: damped , localOrtho , localVQueueFull , prnt , show integer ( ip ) :: i , localOrthoCount , localOrthoLimit , localPointer , localVecs , & pcount , pfreq real ( dp ) :: alpha , alphabar , alphahat , & beta , betaacute , betacheck , betad , betadd , betahat , & normb , c , cbar , chat , ctildeold , ctol , & d , maxrbar , minrbar , normA2 , & rho , rhobar , rhobarold , rhodold , rhoold , rhotemp , & rhotildeold , rtol , s , sbar , shat , stildeold , & t1 , taud , tautildeold , test1 , test2 , test3 , & thetabar , thetanew , thetatilde , thetatildeold , & zeta , zetabar , zetaold ! Local constants real ( dp ), parameter :: zero = 0.0_dp , one = 1.0_dp character ( len =* ), parameter :: enter = ' Enter LSMR. ' character ( len =* ), parameter :: exitt = ' Exit LSMR. ' character ( len =* ), parameter :: msg ( 0 : 7 ) = & ( / 'The exact solution is x = 0 ' , & 'Ax - b is small enough, given atol, btol ' , & 'The least-squares solution is good enough, given atol' , & 'The estimate of cond(Abar) has exceeded conlim ' , & 'Ax - b is small enough for this machine ' , & 'The LS solution is good enough for this machine ' , & 'Cond(Abar) seems to be too large for this machine ' , & 'The iteration limit has been reached ' / ) !------------------------------------------------------------------- ! Initialize. localVecs = min ( localSize , m , n ) show = nout > 0 if ( show ) then write ( nout , 1000 ) enter , m , n , damp , atol , conlim , btol , itnlim , localVecs end if pfreq = 20 ! print frequency (for repeating the heading) pcount = 0 ! print counter damped = damp > zero ! !------------------------------------------------------------------- ! Set up the first vectors u and v for the bidiagonalization. ! These satisfy beta*u = b, alpha*v = A(transpose)*u. !------------------------------------------------------------------- u ( 1 : m ) = b ( 1 : m ) v ( 1 : n ) = zero x ( 1 : n ) = zero alpha = zero beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) call Aprod2 ( m , n , v , u ) ! v = A'*u alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) end if if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) w = v end if normAr = alpha * beta if ( normAr == zero ) go to 800 ! Initialization for local reorthogonalization. localOrtho = . false . if ( localVecs > 0 ) then localPointer = 1 localOrtho = . true . localVQueueFull = . false . localV (:, 1 ) = v end if ! Initialize variables for 1st iteration. itn = 0 zetabar = alpha * beta alphabar = alpha rho = 1 rhobar = 1 cbar = 1 sbar = 0 h = v hbar ( 1 : n ) = zero x ( 1 : n ) = zero ! Initialize variables for estimation of ||r||. betadd = beta betad = 0 rhodold = 1 tautildeold = 0 thetatilde = 0 zeta = 0 d = 0 ! Initialize variables for estimation of ||A|| and cond(A). normA2 = alpha ** 2 maxrbar = 0_dp minrbar = 1e+100_dp ! Items for use in stopping rules. normb = beta istop = 0 ctol = zero if ( conlim > zero ) ctol = one / conlim normr = beta ! Exit if b=0 or A'b = 0. normAr = alpha * beta if ( normAr == 0 ) then if ( show ) then write ( nout , '(a)' ) msg ( 1 ) end if return end if ! Heading for iteration log. if ( show ) then if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if test1 = one test2 = alpha / beta write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 end if !=================================================================== ! Main iteration loop. !=================================================================== do itn = itn + 1 !---------------------------------------------------------------- ! Perform the next step of the bidiagonalization to obtain the ! next beta, u, alpha, v. These satisfy ! beta*u = A*v - alpha*u, ! alpha*v = A'*u - beta*v. !---------------------------------------------------------------- u = ( - alpha ) * u ! call dscal (m,(- alpha), u, 1) call Aprod1 ( m , n , v , u ) ! u = u + A*v beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) if ( localOrtho ) then ! Store v into the circular buffer localV. call localVEnqueue ! Store old v for local reorthog'n of new v. end if v = ( - beta ) * v ! call dscal (n, (- beta), v, 1) call Aprod2 ( m , n , v , u ) ! v = v + A'*u if ( localOrtho ) then ! Perform local reorthogonalization of V. call localVOrtho ! Local-reorthogonalization of new v. end if alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) end if end if ! At this point, beta = beta_{k+1}, alpha = alpha_{k+1}. !---------------------------------------------------------------- ! Construct rotation Qhat_{k,2k+1}. alphahat = d2norm ( alphabar , damp ) chat = alphabar / alphahat shat = damp / alphahat ! Use a plane rotation (Q_i) to turn B_i to R_i. rhoold = rho rho = d2norm ( alphahat , beta ) c = alphahat / rho s = beta / rho thetanew = s * alpha alphabar = c * alpha ! Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar. rhobarold = rhobar zetaold = zeta thetabar = sbar * rho rhotemp = cbar * rho rhobar = d2norm ( cbar * rho , thetanew ) cbar = cbar * rho / rhobar sbar = thetanew / rhobar zeta = cbar * zetabar zetabar = - sbar * zetabar ! Update h, h_hat, x. hbar = h - ( thetabar * rho / ( rhoold * rhobarold )) * hbar x = x + ( zeta / ( rho * rhobar )) * hbar h = v - ( thetanew / rho ) * h ! Estimate ||r||. ! Apply rotation Qhat_{k,2k+1}. betaacute = chat * betadd betacheck = - shat * betadd ! Apply rotation Q_{k,k+1}. betahat = c * betaacute betadd = - s * betaacute ! Apply rotation Qtilde_{k-1}. ! betad = betad_{k-1} here. thetatildeold = thetatilde rhotildeold = d2norm ( rhodold , thetabar ) ctildeold = rhodold / rhotildeold stildeold = thetabar / rhotildeold thetatilde = stildeold * rhobar rhodold = ctildeold * rhobar betad = - stildeold * betad + ctildeold * betahat ! betad = betad_k here. ! rhodold = rhod_k here. tautildeold = ( zetaold - thetatildeold * tautildeold ) / rhotildeold taud = ( zeta - thetatilde * tautildeold ) / rhodold d = d + betacheck ** 2 normr = sqrt ( d + ( betad - taud ) ** 2 + betadd ** 2 ) ! Estimate ||A||. normA2 = normA2 + beta ** 2 normA = sqrt ( normA2 ) normA2 = normA2 + alpha ** 2 ! Estimate cond(A). maxrbar = max ( maxrbar , rhobarold ) if ( itn > 1 ) then minrbar = min ( minrbar , rhobarold ) end if condA = max ( maxrbar , rhotemp ) / min ( minrbar , rhotemp ) !---------------------------------------------------------------- ! Test for convergence. !---------------------------------------------------------------- ! Compute norms for convergence testing. normAr = abs ( zetabar ) normx = sqrt ( dot_product ( x , x ) ) ! dnrm2(n, x, 1) ! Now use these norms to estimate certain other quantities, ! some of which will be small near a solution. test1 = normr / normb test2 = normAr / ( normA * normr ) test3 = one / condA t1 = test1 / ( one + normA * normx / normb ) rtol = btol + atol * normA * normx / normb ! The following tests guard against extremely small values of ! atol, btol or ctol. (The user may have set any or all of ! the parameters atol, btol, conlim to 0.) ! The effect is equivalent to the normAl tests using ! atol = eps, btol = eps, conlim = 1/eps. if ( itn >= itnlim ) istop = 7 if ( one + test3 <= one ) istop = 6 if ( one + test2 <= one ) istop = 5 if ( one + t1 <= one ) istop = 4 ! Allow for tolerances set by the user. if ( test3 <= ctol ) istop = 3 if ( test2 <= atol ) istop = 2 if ( test1 <= rtol ) istop = 1 !---------------------------------------------------------------- ! See if it is time to print something. !---------------------------------------------------------------- prnt = . false . if ( show ) then if ( n <= 40 ) prnt = . true . if ( itn <= 10 ) prnt = . true . if ( itn >= itnlim - 10 ) prnt = . true . if ( mod ( itn , 10 ) == 0 ) prnt = . true . if ( test3 <= 1.1 * ctol ) prnt = . true . if ( test2 <= 1.1 * atol ) prnt = . true . if ( test1 <= 1.1 * rtol ) prnt = . true . if ( istop /= 0 ) prnt = . true . if ( prnt ) then ! Print a line for this iteration if ( pcount >= pfreq ) then ! Print a heading first pcount = 0 if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if end if pcount = pcount + 1 write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 , normA , condA end if end if if ( istop /= 0 ) exit end do !=================================================================== ! End of iteration loop. !=================================================================== ! Come here if normAr = 0, or if normal exit. 800 if ( show ) then ! Print the stopping condition. write ( nout , 2000 ) & exitt , istop , itn , & exitt , normA , condA , & exitt , normb , normx , & exitt , normr , normAr write ( nout , 3000 ) & exitt , msg ( istop ) end if return 1000 format ( // a , ' Least-squares solution of Ax = b' & / ' The matrix A has' , i7 , ' rows and' , i7 , ' columns' & / ' damp =' , es22 . 14 & / ' atol =' , es10 . 2 , 15 x , 'conlim =' , es10 . 2 & / ' btol =' , es10 . 2 , 15 x , 'itnlim =' , i10 & / ' localSize (no. of vectors for local reorthogonalization) =' , i7 ) 1200 format ( / \" Itn x(1) norm r A'r \" , & ' Compatible LS norm A cond A' ) 1300 format ( / \" Itn x(1) norm rbar Abar'rbar\" , & ' Compatible LS norm Abar cond Abar' ) 1500 format ( i6 , 2 es17 . 9 , 5 es10 . 2 ) 2000 format ( / a , 5 x , 'istop =' , i2 , 15 x , 'itn =' , i8 & / a , 5 x , 'normA =' , es12 . 5 , 5 x , 'condA =' , es12 . 5 & / a , 5 x , 'normb =' , es12 . 5 , 5 x , 'normx =' , es12 . 5 & / a , 5 x , 'normr =' , es12 . 5 , 5 x , 'normAr =' , es12 . 5 ) 3000 format ( a , 5 x , a ) contains function d2norm ( a , b ) real ( dp ) :: d2norm real ( dp ), intent ( in ) :: a , b !------------------------------------------------------------------- ! d2norm returns sqrt( a**2 + b**2 ) ! with precautions to avoid overflow. ! ! 21 Mar 1990: First version. ! 17 Sep 2007: Fortran 90 version. ! 24 Oct 2007: User real(dp) instead of compiler option -r8. !------------------------------------------------------------------- intrinsic :: abs , sqrt real ( dp ) :: scale real ( dp ), parameter :: zero = 0.0_dp scale = abs ( a ) + abs ( b ) if ( scale == zero ) then d2norm = zero else d2norm = scale * sqrt (( a / scale ) ** 2 + ( b / scale ) ** 2 ) end if end function d2norm !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVEnqueue ! Store v into the circular buffer localV. if ( localPointer < localVecs ) then localPointer = localPointer + 1 else localPointer = 1 localVQueueFull = . true . end if localV (:, localPointer ) = v end subroutine localVEnqueue !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVOrtho ! Perform local reorthogonalization of current v. real ( dp ) :: d if ( localVQueueFull ) then localOrthoLimit = localVecs else localOrthoLimit = localPointer end if do localOrthoCount = 1 , localOrthoLimit d = dot_product ( v , localV (:, localOrthoCount )) v = v - d * localV (:, localOrthoCount ) end do end subroutine localVOrtho end subroutine lsmr","tags":"","loc":"proc/lsmr.html"},{"title":"lsmr_ez – LSMR","text":"public subroutine lsmr_ez(m, n, irow, icol, a, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) Easy interface to lsmr . Instead of specifying the Aprod1 , Aprod2 functions,\n the sparsity pattern ( irow , icol ) and nonzero elemenets\n of a are input. History JW : 1/24/2024 : created. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n integer, intent(in), dimension(:) :: irow row indices of nonzero elements of A integer, intent(in), dimension(:) :: icol column indices of nonzero elements of A real(kind=dp), intent(in), dimension(:) :: a nonzero elements of A real(kind=dp), intent(in) :: b (m) real(kind=dp), intent(in) :: damp real(kind=dp), intent(in) :: atol real(kind=dp), intent(in) :: btol real(kind=dp), intent(in) :: conlim integer(kind=ip), intent(in) :: itnlim integer(kind=ip), intent(in) :: localSize integer(kind=ip), intent(in) :: nout real(kind=dp), intent(out) :: x (n) integer(kind=ip), intent(out) :: istop integer(kind=ip), intent(out) :: itn real(kind=dp), intent(out) :: normA real(kind=dp), intent(out) :: condA real(kind=dp), intent(out) :: normr real(kind=dp), intent(out) :: normAr real(kind=dp), intent(out) :: normx Calls proc~~lsmr_ez~~CallsGraph proc~lsmr_ez lsmrModule::lsmr_ez proc~lsmr lsmrModule::lsmr proc~lsmr_ez->proc~lsmr Help Graph Key Nodes of different colours represent the following: Graph Key Subroutine Subroutine Function Function Interface Interface Type Bound Procedure Type Bound Procedure Unknown Procedure Type Unknown Procedure Type Program Program This Page's Entity This Page's Entity Solid arrows point from a procedure to one which it calls. Dashed \narrows point from an interface to procedures which implement that interface.\nThis could include the module procedures in a generic interface or the\nimplementation in a submodule of an interface in a parent module. Source Code subroutine lsmr_ez ( m , n , irow , icol , a , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer , dimension (:), intent ( in ) :: irow !! row indices of nonzero elements of `A` integer , dimension (:), intent ( in ) :: icol !! column indices of nonzero elements of `A` real ( dp ), dimension (:), intent ( in ) :: a !! nonzero elements of `A` integer ( ip ), intent ( in ) :: m , n , itnlim , localSize , nout integer ( ip ), intent ( out ) :: istop , itn real ( dp ), intent ( in ) :: b ( m ) real ( dp ), intent ( out ) :: x ( n ) real ( dp ), intent ( in ) :: atol , btol , conlim , damp real ( dp ), intent ( out ) :: normA , condA , normr , normAr , normx integer ( ip ) :: num_nonzero_elements !! number of nonzero elements in the matrix real ( dp ), dimension (:), allocatable :: Ax , Aty !! temp arrays if ( size ( irow ) == size ( icol ) . and . size ( irow ) == size ( a )) then num_nonzero_elements = size ( irow ) allocate ( Ax ( m )) allocate ( Aty ( n )) call lsmr ( m , n , Aprod1_ez , Aprod2_ez , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) else error stop 'inconsistent sizes of input arrays irow, icol, a' end if contains ! see code from LSQR subroutine Aprod1_ez ( m , n , x , y ) !! y := y + A*x integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index ! A*x: Ax = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Ax ( r ) = Ax ( r ) + a ( i ) * x ( c ) end do y = y + Ax end subroutine Aprod1_ez subroutine Aprod2_ez ( m , n , x , y ) !! x := x + A'*y integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index Aty = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Aty ( c ) = Aty ( c ) + a ( i ) * y ( r ) end do x = x + Aty end subroutine Aprod2_ez end subroutine lsmr_ez","tags":"","loc":"proc/lsmr_ez.html"},{"title":"ddot – LSMR","text":"public function ddot(n, dx, incx, dy, incy) Dot product of two vectors. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vectors. real(kind=wp), intent(in) :: dx (*) the first vector integer(kind=ip), intent(in) :: incx the increment between successive entries in DX. real(kind=wp), intent(in) :: dy (*) the second vector integer(kind=ip), intent(in) :: incy the increment between successive entries in DY. Return Value real(kind=wp) Source Code function ddot ( n , dx , incx , dy , incy ) implicit none integer , intent ( in ) :: n , incx , incy double precision , intent ( in ) :: dx ( * ), dy ( * ) double precision :: ddot end function ddot","tags":"","loc":"proc/ddot.html"},{"title":"dnrm2 – LSMR","text":"public function dnrm2(n, x, incx) The euclidean norm of a vector sqrt ( X' * X ) . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n real(kind=wp), intent(in) :: x (*) integer(kind=ip), intent(in) :: incx Return Value real(kind=wp) Source Code function dnrm2 ( n , dx , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: dx ( * ) double precision :: dnrm2 end function dnrm2","tags":"","loc":"proc/dnrm2.html"},{"title":"dcopy – LSMR","text":"public subroutine dcopy(n, dx, incx, dy, incy) Copies a vector X to a vector Y. Arguments Type Intent Optional Attributes Name integer(kind=ip) :: n real(kind=wp) :: dx (*) integer(kind=ip) :: incx real(kind=wp) :: dy (*) integer(kind=ip) :: incy Source Code subroutine dcopy ( n , dx , incx , dy , incy ) implicit none double precision dx ( * ), dy ( * ) integer i , incx , incy , ix , iy , m , n end subroutine dcopy","tags":"","loc":"proc/dcopy.html"},{"title":"dscal – LSMR","text":"public subroutine dscal(n, sa, x, incx) Scales a vector by a constant. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vector. real(kind=wp), intent(in) :: sa the multiplier. real(kind=wp), intent(inout) :: x (*) the vector to be scaled. integer(kind=ip), intent(in) :: incx the increment between successive entries of X. Source Code subroutine dscal ( n , sa , x , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: sa double precision , intent ( inout ) :: x ( * ) end subroutine dscal","tags":"","loc":"proc/dscal.html"},{"title":"lsmrDataModule – LSMR","text":"Defines integer(ip) and real(wp) kind parameters. Notes JW: 1/24/2024 : The original version of this file was rewritten. Uses iso_fortran_env module~~lsmrdatamodule~~UsesGraph module~lsmrdatamodule lsmrDataModule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Used by module~~lsmrdatamodule~~UsedByGraph module~lsmrdatamodule lsmrDataModule module~lsmrblas lsmrblas module~lsmrblas->module~lsmrdatamodule module~lsmrmodule lsmrModule module~lsmrmodule->module~lsmrdatamodule Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Variables Type Visibility Attributes Name Initial integer, public, parameter :: lsmr_wp = real64 real kind used by this module [8 bytes] integer, public, parameter :: lsmr_ip = int32 integer kind used by this module [4 bytes]","tags":"","loc":"module/lsmrdatamodule.html"},{"title":"lsmrModule – LSMR","text":"LSMR solves Ax = b or min ||Ax - b|| with or without damping,\n using the iterative algorithm of David Fong and Michael Saunders. Authors David Fong & Michael Saunders, Systems Optimization Laboratory (SOL) See also http://www.stanford.edu/group/SOL/software/lsmr.html History 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m. 07 Sep 2010: Local reorthogonalization now works (localSize > 0). 28 Jan 2014: In lsmrDataModule.f90:\n ip added for integer(ip) declarations.\n dnrm2 and dscal coded directly\n (no longer use lsmrblasInterface.f90 or lsmrblas.f90). Uses lsmrDataModule module~~lsmrmodule~~UsesGraph module~lsmrmodule lsmrModule module~lsmrdatamodule lsmrDataModule module~lsmrmodule->module~lsmrdatamodule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Abstract Interfaces abstract interface private subroutine Aprod1_f(m, n, x, y) y := y + A*x Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(in) :: x (n) real(kind=dp), intent(inout) :: y (m) abstract interface private subroutine Aprod2_f(m, n, x, y) x := x + A'*y Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(inout) :: x (n) real(kind=dp), intent(in) :: y (m) Subroutines public subroutine lsmr (m, n, Aprod1, Aprod2, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) LSMR finds a solution x to the following problems: Read more… Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m the number of rows in A. integer(kind=ip), intent(in) :: n the number of columns in A. procedure( Aprod1_f ) :: Aprod1 See above. procedure( Aprod2_f ) :: Aprod2 See above. real(kind=dp), intent(in) :: b (m) The rhs vector b . real(kind=dp), intent(in) :: damp The damping parameter for problem 3 above.\n(damp should be 0.0 for problems 1 and 2.)\nIf the system A*x = b is incompatible, values\nof damp in the range 0 to sqrt(eps)*norm(A) will probably have a negligible effect.\nLarger values of damp will tend to decrease\nthe norm of x and reduce the number of\niterations required by LSMR. Read more… real(kind=dp), intent(in) :: atol An estimate of the relative error in the data\ndefining the matrix A. For example, if A is\naccurate to about 6 digits, set atol = 1.0e-6. real(kind=dp), intent(in) :: btol An estimate of the relative error in the data\ndefining the rhs b. For example, if b is\naccurate to about 6 digits, set btol = 1.0e-6. real(kind=dp), intent(in) :: conlim An upper limit on cond(Abar), the apparent\ncondition number of the matrix Abar.\nIterations will be terminated if a computed\nestimate of cond(Abar) exceeds conlim.\nThis is intended to prevent certain small or\nzero singular values of A or Abar from\ncoming into effect and causing unwanted growth\nin the computed solution. Read more… integer(kind=ip), intent(in) :: itnlim An upper limit on the number of iterations.\nSuggested value: Read more… integer(kind=ip), intent(in) :: localSize No. of vectors for local reorthogonalization: Read more… integer(kind=ip), intent(in) :: nout File number for printed output. If positive,\na summary will be printed on file nout. real(kind=dp), intent(out) :: x (n) Returns the computed solution x . integer(kind=ip), intent(out) :: istop An integer giving the reason for termination: Read more… integer(kind=ip), intent(out) :: itn The number of iterations performed. real(kind=dp), intent(out) :: normA An estimate of the Frobenius norm of Abar.\nThis is the square-root of the sum of squares\nof the elements of Abar.\nIf damp is small and the columns of A\nhave all been scaled to have length 1.0,\nnormA should increase to roughly sqrt(n).\nA radically different value for normA may\nindicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: condA An estimate of cond(Abar), the condition\nnumber of Abar. A very high value of condA\nmay again indicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: normr An estimate of the final value of norm(rbar),\nthe function being minimized (see notation\nabove). This will be small if A*x = b has\na solution. real(kind=dp), intent(out) :: normAr An estimate of the final value of norm( Abar'*rbar ) , the norm of\nthe residual for the normal equations.\nThis should be small in all cases. (normAr\nwill often be smaller than the true value\ncomputed from the output vector x.) real(kind=dp), intent(out) :: normx An estimate of norm(x) for the final solution x. public subroutine lsmr_ez (m, n, irow, icol, a, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) Easy interface to lsmr . Read more… Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n integer, intent(in), dimension(:) :: irow row indices of nonzero elements of A integer, intent(in), dimension(:) :: icol column indices of nonzero elements of A real(kind=dp), intent(in), dimension(:) :: a nonzero elements of A real(kind=dp), intent(in) :: b (m) real(kind=dp), intent(in) :: damp real(kind=dp), intent(in) :: atol real(kind=dp), intent(in) :: btol real(kind=dp), intent(in) :: conlim integer(kind=ip), intent(in) :: itnlim integer(kind=ip), intent(in) :: localSize integer(kind=ip), intent(in) :: nout real(kind=dp), intent(out) :: x (n) integer(kind=ip), intent(out) :: istop integer(kind=ip), intent(out) :: itn real(kind=dp), intent(out) :: normA real(kind=dp), intent(out) :: condA real(kind=dp), intent(out) :: normr real(kind=dp), intent(out) :: normAr real(kind=dp), intent(out) :: normx","tags":"","loc":"module/lsmrmodule.html"},{"title":"lsmrblas – LSMR","text":"This file contains the following BLAS routines dcopy , ddot , dnrm2 , dscal required by subroutines lsmr and acheck . References Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,\n LINPACK User's Guide,\n SIAM, 1979,\n ISBN13: 978-0-898711-72-1,\n LC: QA214.L56. Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,\n Algorithm 539,\n Basic Linear Algebra Subprograms for Fortran Usage,\n ACM Transactions on Mathematical Software,\n Volume 5, Number 3, September 1979, pages 308-323. Uses lsmrDataModule module~~lsmrblas~~UsesGraph module~lsmrblas lsmrblas module~lsmrdatamodule lsmrDataModule module~lsmrblas->module~lsmrdatamodule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Functions public function ddot (n, dx, incx, dy, incy) Dot product of two vectors. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vectors. real(kind=wp), intent(in) :: dx (*) the first vector integer(kind=ip), intent(in) :: incx the increment between successive entries in DX. real(kind=wp), intent(in) :: dy (*) the second vector integer(kind=ip), intent(in) :: incy the increment between successive entries in DY. Return Value real(kind=wp) public function dnrm2 (n, x, incx) The euclidean norm of a vector sqrt ( X' * X ) . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n real(kind=wp), intent(in) :: x (*) integer(kind=ip), intent(in) :: incx Return Value real(kind=wp) Subroutines public subroutine dcopy (n, dx, incx, dy, incy) Copies a vector X to a vector Y. Arguments Type Intent Optional Attributes Name integer(kind=ip) :: n real(kind=wp) :: dx (*) integer(kind=ip) :: incx real(kind=wp) :: dy (*) integer(kind=ip) :: incy public subroutine dscal (n, sa, x, incx) Scales a vector by a constant. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vector. real(kind=wp), intent(in) :: sa the multiplier. real(kind=wp), intent(inout) :: x (*) the vector to be scaled. integer(kind=ip), intent(in) :: incx the increment between successive entries of X.","tags":"","loc":"module/lsmrblas.html"},{"title":"lsmrDataModule.F90 – LSMR","text":"Files dependent on this one sourcefile~~lsmrdatamodule.f90~~AfferentGraph sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrblas.f90 lsmrblas.F90 sourcefile~lsmrblas.f90->sourcefile~lsmrdatamodule.f90 sourcefile~lsmrmodule.f90 lsmrModule.f90 sourcefile~lsmrmodule.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! Defines `integer(ip)` and `real(wp)` kind parameters. ! !### Notes ! * JW: 1/24/2024 : The original version of this file was rewritten. module lsmrDataModule use iso_fortran_env implicit none private #ifdef REAL32 integer , parameter , public :: lsmr_wp = real32 !! real kind used by this module [4 bytes] #elif REAL64 integer , parameter , public :: lsmr_wp = real64 !! real kind used by this module [8 bytes] #elif REAL128 integer , parameter , public :: lsmr_wp = real128 !! real kind used by this module [16 bytes] #else integer , parameter , public :: lsmr_wp = real64 !! real kind used by this module [8 bytes] #endif integer , parameter , public :: lsmr_ip = int32 !! integer kind used by this module [4 bytes] end module lsmrDataModule","tags":"","loc":"sourcefile/lsmrdatamodule.f90.html"},{"title":"lsmrModule.f90 – LSMR","text":"This file depends on sourcefile~~lsmrmodule.f90~~EfferentGraph sourcefile~lsmrmodule.f90 lsmrModule.f90 sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrmodule.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! LSMR solves `Ax = b` or `min ||Ax - b||` with or without damping, ! using the iterative algorithm of David Fong and Michael Saunders. ! !### Authors ! * David Fong & Michael Saunders, Systems Optimization Laboratory (SOL) ! !### See also ! * http://www.stanford.edu/group/SOL/software/lsmr.html ! !### History ! * 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m. ! * 07 Sep 2010: Local reorthogonalization now works (localSize > 0). ! * 28 Jan 2014: In lsmrDataModule.f90: ! ip added for integer(ip) declarations. ! dnrm2 and dscal coded directly ! (no longer use lsmrblasInterface.f90 or lsmrblas.f90). module lsmrModule use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none private public :: lsmr public :: lsmr_ez abstract interface subroutine Aprod1_f ( m , n , x , y ) !! y := y + A*x use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) end subroutine Aprod1_f subroutine Aprod2_f ( m , n , x , y ) !! x := x + A'*y use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) end subroutine Aprod2_f end interface contains !***************************************************************************** !> ! LSMR finds a solution x to the following problems: ! ! 1. Unsymmetric equations: Solve `A*x = b` ! 2. Linear least squares: Solve `A*x = b` in the least-squares sense ! 3. Damped least squares: Solve !``` ! ( A )*x = ( b ) ! ( damp*I ) ( 0 ) !``` ! in the least-squares sense ! ! where A is a matrix with m rows and n columns, b is an m-vector, ! and damp is a scalar. (All quantities are real.) ! The matrix A is treated as a linear operator. It is accessed ! by means of subroutine calls with the following purpose: ! ! * `call Aprod1(m,n,x,y)` must compute `y = y + A*x ` without altering `x`. ! * `call Aprod2(m,n,x,y)` must compute `x = x + A'*y` without altering `y`. ! ! LSMR uses an iterative method to approximate the solution. ! The number of iterations required to reach a certain accuracy ! depends strongly on the scaling of the problem. Poor scaling of ! the rows or columns of A should therefore be avoided where ! possible. ! ! For example, in problem 1 the solution is unaltered by ! row-scaling. If a row of A is very small or large compared to ! the other rows of A, the corresponding row of ( A b ) should be ! scaled up or down. ! ! In problems 1 and 2, the solution x is easily recovered ! following column-scaling. Unless better information is known, ! the nonzero columns of A should be scaled so that they all have ! the same Euclidean norm (e.g., 1.0). ! ! In problem 3, there is no freedom to re-scale if damp is ! nonzero. However, the value of damp should be assigned only ! after attention has been paid to the scaling of A. ! ! The parameter damp is intended to help regularize ! ill-conditioned systems, by preventing the true solution from ! being very large. Another aid to regularization is provided by ! the parameter condA, which may be used to terminate iterations ! before the computed solution becomes very large. ! ! Note that x is not an input parameter. ! If some initial estimate x0 is known and if damp = 0, ! one could proceed as follows: ! ! 1. Compute a residual vector `r0 = b - A*x0`. ! 2. Use LSMR to solve the system `A*dx = r0`. ! 3. Add the correction `dx` to obtain a final solution `x = x0 + dx`. ! ! This requires that x0 be available before and after the call ! to LSMR. To judge the benefits, suppose LSMR takes k1 iterations ! to solve `A*x = b` and k2 iterations to solve `A*dx = r0`. ! If x0 is \"good\", norm(r0) will be smaller than norm(b). ! If the same stopping tolerances atol and btol are used for each ! system, k1 and k2 will be similar, but the final solution `x0 + dx` ! should be more accurate. The only way to reduce the total work ! is to use a larger stopping tolerance for the second system. ! If some value btol is suitable for `A*x = b`, the larger value ! `btol*norm(b)/norm(r0)` should be suitable for `A*dx = r0`. ! ! Preconditioning is another way to reduce the number of iterations. ! If it is possible to solve a related system `M*x = b` efficiently, ! where M approximates A in some helpful way ! (e.g. M - A has low rank or its elements are small relative to ! those of A), LSMR may converge more rapidly on the system ! `A*M(inverse)*z = b`, ! after which `x` can be recovered by solving `M*x = z`. ! ! NOTE: If A is symmetric, LSMR should not be used! ! Alternatives are the symmetric conjugate-gradient method (CG) ! and/or SYMMLQ. ! SYMMLQ is an implementation of symmetric CG that applies to ! any symmetric A and will converge more rapidly than LSMR. ! If A is positive definite, there are other implementations of ! symmetric CG that require slightly less work per iteration ! than SYMMLQ (but will take the same number of iterations). ! ! Notation ! -------- ! The following quantities are used in discussing the subroutine ! parameters: !``` ! Abar = ( A ), bbar = (b) ! (damp*I) (0) ! ! r = b - A*x, rbar = bbar - Abar*x ! ! normr = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) ! = norm( rbar ) ! ! eps = the relative precision of floating-point arithmetic. ! On most machines, eps is about 1.0e-7 and 1.0e-16 ! in single and double precision respectively. ! We expect eps to be about 1e-16 always. !``` ! LSMR minimizes the function `normr` with respect to `x`. ! !### Precision ! ! The number of iterations required by LSMR will decrease ! if the computation is performed in higher precision. ! !### Reference ! * http://www.stanford.edu/group/SOL/software/lsmr.html ! !### LSMR development: ! * 21 Sep 2007: Fortran 90 version of LSQR implemented. ! Aprod1, Aprod2 implemented via f90 interface. ! * 17 Jul 2010: LSMR derived from LSQR and lsmr.m. ! * 07 Sep 2010: Local reorthogonalization now working. ! * 02 May 2014: With damp>0, istop=2 was incorrectly set to istop=3 ! (so incorrect stopping message was printed). Fixed. ! !@note Any or all of `atol`, `btol`, `conlim` may be set to zero. ! The effect will be the same as the values `eps`, `eps`, `1/eps`. subroutine lsmr ( m , n , Aprod1 , Aprod2 , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer ( ip ), intent ( in ) :: m !! the number of rows in A. integer ( ip ), intent ( in ) :: n !! the number of columns in A. integer ( ip ), intent ( in ) :: itnlim !! An upper limit on the number of iterations. !! Suggested value: !! !! * `itnlim = n/2` for well-conditioned systems !! with clustered singular values, !! * `itnlim = 4*n` otherwise. integer ( ip ), intent ( in ) :: localSize !! No. of vectors for local reorthogonalization: !! !! * 0 No reorthogonalization is performed. !! * >0 This many n-vectors \"v\" (the most recent ones) !! are saved for reorthogonalizing the next v. !! !! localSize need not be more than min(m,n). !! At most min(m,n) vectors will be allocated. integer ( ip ), intent ( in ) :: nout !! File number for printed output. If positive, !! a summary will be printed on file nout. integer ( ip ), intent ( out ) :: istop !! An integer giving the reason for termination: !! !! * 0 x = 0 is the exact solution. !! No iterations were performed. !! * 1 The equations `A*x = b` are probably compatible. !! `Norm(A*x - b)` is sufficiently small, given the !! values of atol and btol. !! * 2 damp is zero. The system `A*x = b` is probably !! not compatible. A least-squares solution has !! been obtained that is sufficiently accurate, !! given the value of atol. !! * 3 damp is nonzero. A damped least-squares !! solution has been obtained that is sufficiently !! accurate, given the value of atol. !! * 4 An estimate of cond(Abar) has exceeded conlim. !! The system `A*x = b` appears to be ill-conditioned, !! or there could be an error in Aprod1 or Aprod2. !! * 5 The iteration limit itnlim was reached. integer ( ip ), intent ( out ) :: itn !! The number of iterations performed. real ( dp ), intent ( in ) :: b ( m ) !! The rhs vector `b`. real ( dp ), intent ( out ) :: x ( n ) !! Returns the computed solution `x`. real ( dp ), intent ( in ) :: atol !! An estimate of the relative error in the data !! defining the matrix A. For example, if A is !! accurate to about 6 digits, set atol = 1.0e-6. real ( dp ), intent ( in ) :: btol !! An estimate of the relative error in the data !! defining the rhs b. For example, if b is !! accurate to about 6 digits, set btol = 1.0e-6. real ( dp ), intent ( in ) :: conlim !! An upper limit on cond(Abar), the apparent !! condition number of the matrix Abar. !! Iterations will be terminated if a computed !! estimate of cond(Abar) exceeds conlim. !! This is intended to prevent certain small or !! zero singular values of A or Abar from !! coming into effect and causing unwanted growth !! in the computed solution. !! !! conlim and damp may be used separately or !! together to regularize ill-conditioned systems. !! !! Normally, conlim should be in the range !! 1000 to 1/eps. !! Suggested value: !! !! * `conlim = 1/(100*eps)` for compatible systems, !! * `conlim = 1/(10*sqrt(eps))` for least squares. real ( dp ), intent ( in ) :: damp !! The damping parameter for problem 3 above. !! (damp should be 0.0 for problems 1 and 2.) !! If the system `A*x = b` is incompatible, values !! of damp in the range 0 to `sqrt(eps)*norm(A)` !! will probably have a negligible effect. !! Larger values of damp will tend to decrease !! the norm of x and reduce the number of !! iterations required by LSMR. !! !! The work per iteration and the storage needed !! by LSMR are the same for all values of damp. real ( dp ), intent ( out ) :: normA !! An estimate of the Frobenius norm of Abar. !! This is the square-root of the sum of squares !! of the elements of Abar. !! If damp is small and the columns of A !! have all been scaled to have length 1.0, !! normA should increase to roughly sqrt(n). !! A radically different value for normA may !! indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: condA !! An estimate of cond(Abar), the condition !! number of Abar. A very high value of condA !! may again indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: normr !! An estimate of the final value of norm(rbar), !! the function being minimized (see notation !! above). This will be small if A*x = b has !! a solution. real ( dp ), intent ( out ) :: normAr !! An estimate of the final value of !! `norm( Abar'*rbar )`, the norm of !! the residual for the normal equations. !! This should be small in all cases. (normAr !! will often be smaller than the true value !! computed from the output vector x.) real ( dp ), intent ( out ) :: normx !! An estimate of norm(x) for the final solution x. procedure ( Aprod1_f ) :: Aprod1 !! See above. procedure ( Aprod2_f ) :: Aprod2 !! See above. ! Local arrays and variables real ( dp ) :: h ( n ), hbar ( n ), u ( m ), v ( n ), w ( n ), localV ( n , min ( localSize , m , n )) logical :: damped , localOrtho , localVQueueFull , prnt , show integer ( ip ) :: i , localOrthoCount , localOrthoLimit , localPointer , localVecs , & pcount , pfreq real ( dp ) :: alpha , alphabar , alphahat , & beta , betaacute , betacheck , betad , betadd , betahat , & normb , c , cbar , chat , ctildeold , ctol , & d , maxrbar , minrbar , normA2 , & rho , rhobar , rhobarold , rhodold , rhoold , rhotemp , & rhotildeold , rtol , s , sbar , shat , stildeold , & t1 , taud , tautildeold , test1 , test2 , test3 , & thetabar , thetanew , thetatilde , thetatildeold , & zeta , zetabar , zetaold ! Local constants real ( dp ), parameter :: zero = 0.0_dp , one = 1.0_dp character ( len =* ), parameter :: enter = ' Enter LSMR. ' character ( len =* ), parameter :: exitt = ' Exit LSMR. ' character ( len =* ), parameter :: msg ( 0 : 7 ) = & ( / 'The exact solution is x = 0 ' , & 'Ax - b is small enough, given atol, btol ' , & 'The least-squares solution is good enough, given atol' , & 'The estimate of cond(Abar) has exceeded conlim ' , & 'Ax - b is small enough for this machine ' , & 'The LS solution is good enough for this machine ' , & 'Cond(Abar) seems to be too large for this machine ' , & 'The iteration limit has been reached ' / ) !------------------------------------------------------------------- ! Initialize. localVecs = min ( localSize , m , n ) show = nout > 0 if ( show ) then write ( nout , 1000 ) enter , m , n , damp , atol , conlim , btol , itnlim , localVecs end if pfreq = 20 ! print frequency (for repeating the heading) pcount = 0 ! print counter damped = damp > zero ! !------------------------------------------------------------------- ! Set up the first vectors u and v for the bidiagonalization. ! These satisfy beta*u = b, alpha*v = A(transpose)*u. !------------------------------------------------------------------- u ( 1 : m ) = b ( 1 : m ) v ( 1 : n ) = zero x ( 1 : n ) = zero alpha = zero beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) call Aprod2 ( m , n , v , u ) ! v = A'*u alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) end if if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) w = v end if normAr = alpha * beta if ( normAr == zero ) go to 800 ! Initialization for local reorthogonalization. localOrtho = . false . if ( localVecs > 0 ) then localPointer = 1 localOrtho = . true . localVQueueFull = . false . localV (:, 1 ) = v end if ! Initialize variables for 1st iteration. itn = 0 zetabar = alpha * beta alphabar = alpha rho = 1 rhobar = 1 cbar = 1 sbar = 0 h = v hbar ( 1 : n ) = zero x ( 1 : n ) = zero ! Initialize variables for estimation of ||r||. betadd = beta betad = 0 rhodold = 1 tautildeold = 0 thetatilde = 0 zeta = 0 d = 0 ! Initialize variables for estimation of ||A|| and cond(A). normA2 = alpha ** 2 maxrbar = 0_dp minrbar = 1e+100_dp ! Items for use in stopping rules. normb = beta istop = 0 ctol = zero if ( conlim > zero ) ctol = one / conlim normr = beta ! Exit if b=0 or A'b = 0. normAr = alpha * beta if ( normAr == 0 ) then if ( show ) then write ( nout , '(a)' ) msg ( 1 ) end if return end if ! Heading for iteration log. if ( show ) then if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if test1 = one test2 = alpha / beta write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 end if !=================================================================== ! Main iteration loop. !=================================================================== do itn = itn + 1 !---------------------------------------------------------------- ! Perform the next step of the bidiagonalization to obtain the ! next beta, u, alpha, v. These satisfy ! beta*u = A*v - alpha*u, ! alpha*v = A'*u - beta*v. !---------------------------------------------------------------- u = ( - alpha ) * u ! call dscal (m,(- alpha), u, 1) call Aprod1 ( m , n , v , u ) ! u = u + A*v beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) if ( localOrtho ) then ! Store v into the circular buffer localV. call localVEnqueue ! Store old v for local reorthog'n of new v. end if v = ( - beta ) * v ! call dscal (n, (- beta), v, 1) call Aprod2 ( m , n , v , u ) ! v = v + A'*u if ( localOrtho ) then ! Perform local reorthogonalization of V. call localVOrtho ! Local-reorthogonalization of new v. end if alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) end if end if ! At this point, beta = beta_{k+1}, alpha = alpha_{k+1}. !---------------------------------------------------------------- ! Construct rotation Qhat_{k,2k+1}. alphahat = d2norm ( alphabar , damp ) chat = alphabar / alphahat shat = damp / alphahat ! Use a plane rotation (Q_i) to turn B_i to R_i. rhoold = rho rho = d2norm ( alphahat , beta ) c = alphahat / rho s = beta / rho thetanew = s * alpha alphabar = c * alpha ! Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar. rhobarold = rhobar zetaold = zeta thetabar = sbar * rho rhotemp = cbar * rho rhobar = d2norm ( cbar * rho , thetanew ) cbar = cbar * rho / rhobar sbar = thetanew / rhobar zeta = cbar * zetabar zetabar = - sbar * zetabar ! Update h, h_hat, x. hbar = h - ( thetabar * rho / ( rhoold * rhobarold )) * hbar x = x + ( zeta / ( rho * rhobar )) * hbar h = v - ( thetanew / rho ) * h ! Estimate ||r||. ! Apply rotation Qhat_{k,2k+1}. betaacute = chat * betadd betacheck = - shat * betadd ! Apply rotation Q_{k,k+1}. betahat = c * betaacute betadd = - s * betaacute ! Apply rotation Qtilde_{k-1}. ! betad = betad_{k-1} here. thetatildeold = thetatilde rhotildeold = d2norm ( rhodold , thetabar ) ctildeold = rhodold / rhotildeold stildeold = thetabar / rhotildeold thetatilde = stildeold * rhobar rhodold = ctildeold * rhobar betad = - stildeold * betad + ctildeold * betahat ! betad = betad_k here. ! rhodold = rhod_k here. tautildeold = ( zetaold - thetatildeold * tautildeold ) / rhotildeold taud = ( zeta - thetatilde * tautildeold ) / rhodold d = d + betacheck ** 2 normr = sqrt ( d + ( betad - taud ) ** 2 + betadd ** 2 ) ! Estimate ||A||. normA2 = normA2 + beta ** 2 normA = sqrt ( normA2 ) normA2 = normA2 + alpha ** 2 ! Estimate cond(A). maxrbar = max ( maxrbar , rhobarold ) if ( itn > 1 ) then minrbar = min ( minrbar , rhobarold ) end if condA = max ( maxrbar , rhotemp ) / min ( minrbar , rhotemp ) !---------------------------------------------------------------- ! Test for convergence. !---------------------------------------------------------------- ! Compute norms for convergence testing. normAr = abs ( zetabar ) normx = sqrt ( dot_product ( x , x ) ) ! dnrm2(n, x, 1) ! Now use these norms to estimate certain other quantities, ! some of which will be small near a solution. test1 = normr / normb test2 = normAr / ( normA * normr ) test3 = one / condA t1 = test1 / ( one + normA * normx / normb ) rtol = btol + atol * normA * normx / normb ! The following tests guard against extremely small values of ! atol, btol or ctol. (The user may have set any or all of ! the parameters atol, btol, conlim to 0.) ! The effect is equivalent to the normAl tests using ! atol = eps, btol = eps, conlim = 1/eps. if ( itn >= itnlim ) istop = 7 if ( one + test3 <= one ) istop = 6 if ( one + test2 <= one ) istop = 5 if ( one + t1 <= one ) istop = 4 ! Allow for tolerances set by the user. if ( test3 <= ctol ) istop = 3 if ( test2 <= atol ) istop = 2 if ( test1 <= rtol ) istop = 1 !---------------------------------------------------------------- ! See if it is time to print something. !---------------------------------------------------------------- prnt = . false . if ( show ) then if ( n <= 40 ) prnt = . true . if ( itn <= 10 ) prnt = . true . if ( itn >= itnlim - 10 ) prnt = . true . if ( mod ( itn , 10 ) == 0 ) prnt = . true . if ( test3 <= 1.1 * ctol ) prnt = . true . if ( test2 <= 1.1 * atol ) prnt = . true . if ( test1 <= 1.1 * rtol ) prnt = . true . if ( istop /= 0 ) prnt = . true . if ( prnt ) then ! Print a line for this iteration if ( pcount >= pfreq ) then ! Print a heading first pcount = 0 if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if end if pcount = pcount + 1 write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 , normA , condA end if end if if ( istop /= 0 ) exit end do !=================================================================== ! End of iteration loop. !=================================================================== ! Come here if normAr = 0, or if normal exit. 800 if ( show ) then ! Print the stopping condition. write ( nout , 2000 ) & exitt , istop , itn , & exitt , normA , condA , & exitt , normb , normx , & exitt , normr , normAr write ( nout , 3000 ) & exitt , msg ( istop ) end if return 1000 format ( // a , ' Least-squares solution of Ax = b' & / ' The matrix A has' , i7 , ' rows and' , i7 , ' columns' & / ' damp =' , es22 . 14 & / ' atol =' , es10 . 2 , 15 x , 'conlim =' , es10 . 2 & / ' btol =' , es10 . 2 , 15 x , 'itnlim =' , i10 & / ' localSize (no. of vectors for local reorthogonalization) =' , i7 ) 1200 format ( / \" Itn x(1) norm r A'r \" , & ' Compatible LS norm A cond A' ) 1300 format ( / \" Itn x(1) norm rbar Abar'rbar\" , & ' Compatible LS norm Abar cond Abar' ) 1500 format ( i6 , 2 es17 . 9 , 5 es10 . 2 ) 2000 format ( / a , 5 x , 'istop =' , i2 , 15 x , 'itn =' , i8 & / a , 5 x , 'normA =' , es12 . 5 , 5 x , 'condA =' , es12 . 5 & / a , 5 x , 'normb =' , es12 . 5 , 5 x , 'normx =' , es12 . 5 & / a , 5 x , 'normr =' , es12 . 5 , 5 x , 'normAr =' , es12 . 5 ) 3000 format ( a , 5 x , a ) contains function d2norm ( a , b ) real ( dp ) :: d2norm real ( dp ), intent ( in ) :: a , b !------------------------------------------------------------------- ! d2norm returns sqrt( a**2 + b**2 ) ! with precautions to avoid overflow. ! ! 21 Mar 1990: First version. ! 17 Sep 2007: Fortran 90 version. ! 24 Oct 2007: User real(dp) instead of compiler option -r8. !------------------------------------------------------------------- intrinsic :: abs , sqrt real ( dp ) :: scale real ( dp ), parameter :: zero = 0.0_dp scale = abs ( a ) + abs ( b ) if ( scale == zero ) then d2norm = zero else d2norm = scale * sqrt (( a / scale ) ** 2 + ( b / scale ) ** 2 ) end if end function d2norm !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVEnqueue ! Store v into the circular buffer localV. if ( localPointer < localVecs ) then localPointer = localPointer + 1 else localPointer = 1 localVQueueFull = . true . end if localV (:, localPointer ) = v end subroutine localVEnqueue !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVOrtho ! Perform local reorthogonalization of current v. real ( dp ) :: d if ( localVQueueFull ) then localOrthoLimit = localVecs else localOrthoLimit = localPointer end if do localOrthoCount = 1 , localOrthoLimit d = dot_product ( v , localV (:, localOrthoCount )) v = v - d * localV (:, localOrthoCount ) end do end subroutine localVOrtho end subroutine lsmr !***************************************************************************** !> ! Easy interface to [[lsmr]]. ! ! Instead of specifying the `Aprod1`, `Aprod2` functions, ! the sparsity pattern (`irow`, `icol`) and nonzero elemenets ! of `a` are input. ! !### History ! * JW : 1/24/2024 : created. subroutine lsmr_ez ( m , n , irow , icol , a , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer , dimension (:), intent ( in ) :: irow !! row indices of nonzero elements of `A` integer , dimension (:), intent ( in ) :: icol !! column indices of nonzero elements of `A` real ( dp ), dimension (:), intent ( in ) :: a !! nonzero elements of `A` integer ( ip ), intent ( in ) :: m , n , itnlim , localSize , nout integer ( ip ), intent ( out ) :: istop , itn real ( dp ), intent ( in ) :: b ( m ) real ( dp ), intent ( out ) :: x ( n ) real ( dp ), intent ( in ) :: atol , btol , conlim , damp real ( dp ), intent ( out ) :: normA , condA , normr , normAr , normx integer ( ip ) :: num_nonzero_elements !! number of nonzero elements in the matrix real ( dp ), dimension (:), allocatable :: Ax , Aty !! temp arrays if ( size ( irow ) == size ( icol ) . and . size ( irow ) == size ( a )) then num_nonzero_elements = size ( irow ) allocate ( Ax ( m )) allocate ( Aty ( n )) call lsmr ( m , n , Aprod1_ez , Aprod2_ez , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) else error stop 'inconsistent sizes of input arrays irow, icol, a' end if contains ! see code from LSQR subroutine Aprod1_ez ( m , n , x , y ) !! y := y + A*x integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index ! A*x: Ax = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Ax ( r ) = Ax ( r ) + a ( i ) * x ( c ) end do y = y + Ax end subroutine Aprod1_ez subroutine Aprod2_ez ( m , n , x , y ) !! x := x + A'*y integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index Aty = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Aty ( c ) = Aty ( c ) + a ( i ) * y ( r ) end do x = x + Aty end subroutine Aprod2_ez end subroutine lsmr_ez !***************************************************************************** end module LSMRmodule","tags":"","loc":"sourcefile/lsmrmodule.f90.html"},{"title":"lsmrblas.F90 – LSMR","text":"This file depends on sourcefile~~lsmrblas.f90~~EfferentGraph sourcefile~lsmrblas.f90 lsmrblas.F90 sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrblas.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! This file contains the following BLAS routines ! [[dcopy]], [[ddot]], [[dnrm2]], [[dscal]] ! required by subroutines [[lsmr]] and [[acheck]]. ! !### References ! ! * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. module lsmrblas use lsmrDataModule , only : ip => lsmr_ip , wp => lsmr_wp implicit none private public :: ddot , dnrm2 , dscal , dcopy #ifdef HAS_BLAS ! get BLAS from an externally-linked library (double precision only) interface function ddot ( n , dx , incx , dy , incy ) implicit none integer , intent ( in ) :: n , incx , incy double precision , intent ( in ) :: dx ( * ), dy ( * ) double precision :: ddot end function ddot function dnrm2 ( n , dx , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: dx ( * ) double precision :: dnrm2 end function dnrm2 subroutine dscal ( n , sa , x , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: sa double precision , intent ( inout ) :: x ( * ) end subroutine dscal subroutine dcopy ( n , dx , incx , dy , incy ) implicit none double precision dx ( * ), dy ( * ) integer i , incx , incy , ix , iy , m , n end subroutine dcopy end interface #else contains !***************************************************************************** !> ! Copies a vector X to a vector Y. subroutine dcopy ( n , dx , incx , dy , incy ) implicit none real ( wp ) dx ( * ), dy ( * ) integer ( ip ) i , incx , incy , ix , iy , m , n if ( n <= 0 ) then return end if if ( incx == 1 . and . incy == 1 ) then m = mod ( n , 7 ) if ( m /= 0 ) then dy ( 1 : m ) = dx ( 1 : m ) end if do i = m + 1 , n , 7 dy ( i ) = dx ( i ) dy ( i + 1 ) = dx ( i + 1 ) dy ( i + 2 ) = dx ( i + 2 ) dy ( i + 3 ) = dx ( i + 3 ) dy ( i + 4 ) = dx ( i + 4 ) dy ( i + 5 ) = dx ( i + 5 ) dy ( i + 6 ) = dx ( i + 6 ) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1 , n dy ( iy ) = dx ( ix ) ix = ix + incx iy = iy + incy end do end if end subroutine dcopy !***************************************************************************** !> ! Dot product of two vectors. real ( wp ) function ddot ( n , dx , incx , dy , incy ) integer ( ip ), intent ( in ) :: n !! the number of entries in the vectors. real ( wp ), intent ( in ) :: dx ( * ) !! the first vector integer ( ip ), intent ( in ) :: incx !! the increment between successive entries in DX. real ( wp ), intent ( in ) :: dy ( * ) !! the second vector integer ( ip ), intent ( in ) :: incy !! the increment between successive entries in DY. integer ( ip ) :: i , ix , iy , m real ( wp ) :: dtemp ddot = 0.0_wp dtemp = 0.0_wp if ( n <= 0 ) return if ( incx /= 1 . or . incy /= 1 ) then ! Code for unequal increments or equal increments ! not equal to 1. if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1 , n dtemp = dtemp + dx ( ix ) * dy ( iy ) ix = ix + incx iy = iy + incy end do else ! Code for both increments equal to 1. m = mod ( n , 5 ) do i = 1 , m dtemp = dtemp + dx ( i ) * dy ( i ) end do do i = m + 1 , n , 5 dtemp = dtemp + dx ( i ) * dy ( i ) + dx ( i + 1 ) * dy ( i + 1 ) + dx ( i + 2 ) * dy ( i + 2 ) & + dx ( i + 3 ) * dy ( i + 3 ) + dx ( i + 4 ) * dy ( i + 4 ) end do end if ddot = dtemp end function ddot !***************************************************************************** !> ! The euclidean norm of a vector `sqrt ( X' * X )`. real ( wp ) function dnrm2 ( n , x , incx ) integer ( ip ), intent ( in ) :: n real ( wp ), intent ( in ) :: x ( * ) integer ( ip ), intent ( in ) :: incx integer ( ip ) :: ix real ( wp ) :: ssq , absxi , norm , scale if ( n < 1 . or . incx < 1 ) then norm = 0.0_wp else if ( n == 1 ) then norm = abs ( x ( 1 ) ) else scale = 0.0_wp ssq = 1.0_wp do ix = 1 , 1 + ( n - 1 ) * incx , incx if ( x ( ix ) /= 0.0_wp ) then absxi = abs ( x ( ix ) ) if ( scale < absxi ) then ssq = 1.0_wp + ssq * ( scale / absxi ) ** 2 scale = absxi else ssq = ssq + ( absxi / scale ) ** 2 end if end if end do norm = scale * sqrt ( ssq ) end if dnrm2 = norm end function dnrm2 !***************************************************************************** !> ! Scales a vector by a constant. subroutine dscal ( n , sa , x , incx ) integer ( ip ), intent ( in ) :: n !! the number of entries in the vector. real ( wp ) , intent ( in ) :: sa !! the multiplier. real ( wp ), intent ( inout ) :: x ( * ) !! the vector to be scaled. integer ( ip ), intent ( in ) :: incx !! the increment between successive entries of X. integer ( ip ) :: i , ix integer ( ip ) :: m if ( n <= 0 ) then return else if ( incx == 1 ) then m = mod ( n , 5 ) x ( 1 : m ) = sa * x ( 1 : m ) do i = m + 1 , n , 5 x ( i ) = sa * x ( i ) x ( i + 1 ) = sa * x ( i + 1 ) x ( i + 2 ) = sa * x ( i + 2 ) x ( i + 3 ) = sa * x ( i + 3 ) x ( i + 4 ) = sa * x ( i + 4 ) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if do i = 1 , n x ( ix ) = sa * x ( ix ) ix = ix + incx end do end if end subroutine dscal #endif end module lsmrblas","tags":"","loc":"sourcefile/lsmrblas.f90.html"}]} \ No newline at end of file +var tipuesearch = {"pages":[{"title":" LSMR ","text":"LSMR LSMR: Sparse Equations and Least Squares. This is a slightly modernized version of the original FSMR (f90 version) by David Fong & Michael Saunders, Systems Optimization Laboratory (SOL), Stanford University, Stanford, CA 94305-4026, USA. Status Compiling A Fortran Package Manager manifest file is included, so that the library and test cases can be compiled with FPM. For example: fpm build --profile release\nfpm test --profile release To use lsmr within your fpm project, add the following to your fpm.toml file: [dependencies] LSMR = { git = \"https://github.com/jacobwilliams/LSMR.git\" } License The original version of LSMR (f90 version) was provided by SOL, Stanford University under the terms of the OSI Common Public License (CPL) or the BSD License . Documentation The latest API documentation can be found here . This was generated from the source code using FORD . See also LSMR: Sparse Equations and Least Squares -- LSMR original code. Sparse Parallel Robust Algorithms Library -- another version Developer Info Jacob Williams","tags":"home","loc":"index.html"},{"title":"Aprod1_f – LSMR","text":"interface private subroutine Aprod1_f(m, n, x, y) Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(in) :: x (n) real(kind=dp), intent(inout) :: y (m) Description y := y + A*x","tags":"","loc":"interface/aprod1_f.html"},{"title":"Aprod2_f – LSMR","text":"interface private subroutine Aprod2_f(m, n, x, y) Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(inout) :: x (n) real(kind=dp), intent(in) :: y (m) Description x := x + A'*y","tags":"","loc":"interface/aprod2_f.html"},{"title":"lsmr – LSMR","text":"public subroutine lsmr(m, n, Aprod1, Aprod2, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) LSMR finds a solution x to the following problems: Unsymmetric equations: Solve A*x = b Linear least squares: Solve A*x = b in the least-squares sense Damped least squares: Solve ( A )*x = ( b )\n ( damp*I ) ( 0 ) in the least-squares sense where A is a matrix with m rows and n columns, b is an m-vector,\n and damp is a scalar. (All quantities are real.)\n The matrix A is treated as a linear operator. It is accessed\n by means of subroutine calls with the following purpose: call Aprod1(m,n,x,y) must compute y = y + A*x without altering x . call Aprod2(m,n,x,y) must compute x = x + A'*y without altering y . LSMR uses an iterative method to approximate the solution.\n The number of iterations required to reach a certain accuracy\n depends strongly on the scaling of the problem. Poor scaling of\n the rows or columns of A should therefore be avoided where\n possible. For example, in problem 1 the solution is unaltered by\n row-scaling. If a row of A is very small or large compared to\n the other rows of A, the corresponding row of ( A b ) should be\n scaled up or down. In problems 1 and 2, the solution x is easily recovered\n following column-scaling. Unless better information is known,\n the nonzero columns of A should be scaled so that they all have\n the same Euclidean norm (e.g., 1.0). In problem 3, there is no freedom to re-scale if damp is\n nonzero. However, the value of damp should be assigned only\n after attention has been paid to the scaling of A. The parameter damp is intended to help regularize\n ill-conditioned systems, by preventing the true solution from\n being very large. Another aid to regularization is provided by\n the parameter condA, which may be used to terminate iterations\n before the computed solution becomes very large. Note that x is not an input parameter.\n If some initial estimate x0 is known and if damp = 0,\n one could proceed as follows: Compute a residual vector r0 = b - A*x0 . Use LSMR to solve the system A*dx = r0 . Add the correction dx to obtain a final solution x = x0 + dx . This requires that x0 be available before and after the call\n to LSMR. To judge the benefits, suppose LSMR takes k1 iterations\n to solve A*x = b and k2 iterations to solve A*dx = r0 .\n If x0 is \"good\", norm(r0) will be smaller than norm(b).\n If the same stopping tolerances atol and btol are used for each\n system, k1 and k2 will be similar, but the final solution x0 + dx should be more accurate. The only way to reduce the total work\n is to use a larger stopping tolerance for the second system.\n If some value btol is suitable for A*x = b , the larger value btol*norm(b)/norm(r0) should be suitable for A*dx = r0 . Preconditioning is another way to reduce the number of iterations.\n If it is possible to solve a related system M*x = b efficiently,\n where M approximates A in some helpful way\n (e.g. M - A has low rank or its elements are small relative to\n those of A), LSMR may converge more rapidly on the system A*M(inverse)*z = b ,\n after which x can be recovered by solving M*x = z . NOTE: If A is symmetric, LSMR should not be used!\n Alternatives are the symmetric conjugate-gradient method (CG)\n and/or SYMMLQ.\n SYMMLQ is an implementation of symmetric CG that applies to\n any symmetric A and will converge more rapidly than LSMR.\n If A is positive definite, there are other implementations of\n symmetric CG that require slightly less work per iteration\n than SYMMLQ (but will take the same number of iterations). Notation The following quantities are used in discussing the subroutine\n parameters: Abar = ( A ), bbar = (b)\n (damp*I) (0)\n\n r = b - A*x, rbar = bbar - Abar*x\n\n normr = sqrt( norm(r)**2 + damp**2 * norm(x)* *2 )\n = norm( rbar )\n\n eps = the relative precision of floating-point arithmetic.\n On most machines, eps is about 1.0e-7 and 1.0e-16\n in single and double precision respectively.\n We expect eps to be about 1e-16 always. LSMR minimizes the function normr with respect to x . Precision The number of iterations required by LSMR will decrease\n if the computation is performed in higher precision. Reference http://www.stanford.edu/group/SOL/software/lsmr.html LSMR development: 21 Sep 2007: Fortran 90 version of LSQR implemented.\n Aprod1, Aprod2 implemented via f90 interface. 17 Jul 2010: LSMR derived from LSQR and lsmr.m. 07 Sep 2010: Local reorthogonalization now working. 02 May 2014: With damp>0, istop=2 was incorrectly set to istop=3\n (so incorrect stopping message was printed). Fixed. Note Any or all of atol , btol , conlim may be set to zero.\n The effect will be the same as the values eps , eps , 1/eps . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m the number of rows in A. integer(kind=ip), intent(in) :: n the number of columns in A. procedure( Aprod1_f ) :: Aprod1 See above. procedure( Aprod2_f ) :: Aprod2 See above. real(kind=dp), intent(in) :: b (m) The rhs vector b . real(kind=dp), intent(in) :: damp The damping parameter for problem 3 above.\n(damp should be 0.0 for problems 1 and 2.)\nIf the system A*x = b is incompatible, values\nof damp in the range 0 to sqrt(eps)*norm(A) will probably have a negligible effect.\nLarger values of damp will tend to decrease\nthe norm of x and reduce the number of\niterations required by LSMR. The work per iteration and the storage needed\nby LSMR are the same for all values of damp. real(kind=dp), intent(in) :: atol An estimate of the relative error in the data\ndefining the matrix A. For example, if A is\naccurate to about 6 digits, set atol = 1.0e-6. real(kind=dp), intent(in) :: btol An estimate of the relative error in the data\ndefining the rhs b. For example, if b is\naccurate to about 6 digits, set btol = 1.0e-6. real(kind=dp), intent(in) :: conlim An upper limit on cond(Abar), the apparent\ncondition number of the matrix Abar.\nIterations will be terminated if a computed\nestimate of cond(Abar) exceeds conlim.\nThis is intended to prevent certain small or\nzero singular values of A or Abar from\ncoming into effect and causing unwanted growth\nin the computed solution. conlim and damp may be used separately or\ntogether to regularize ill-conditioned systems. Normally, conlim should be in the range\n1000 to 1/eps.\nSuggested value: conlim = 1/(100*eps) for compatible systems, conlim = 1/(10*sqrt(eps)) for least squares. integer(kind=ip), intent(in) :: itnlim An upper limit on the number of iterations.\nSuggested value: itnlim = n/2 for well-conditioned systems\n with clustered singular values, itnlim = 4*n otherwise. integer(kind=ip), intent(in) :: localSize No. of vectors for local reorthogonalization: 0 No reorthogonalization is performed. 0 This many n-vectors \"v\" (the most recent ones)\n are saved for reorthogonalizing the next v. localSize need not be more than min(m,n).\nAt most min(m,n) vectors will be allocated. integer(kind=ip), intent(in) :: nout File number for printed output. If positive,\na summary will be printed on file nout. real(kind=dp), intent(out) :: x (n) Returns the computed solution x . integer(kind=ip), intent(out) :: istop An integer giving the reason for termination: 0 x = 0 is the exact solution.\n No iterations were performed. 1 The equations A*x = b are probably compatible. Norm(A*x - b) is sufficiently small, given the\n values of atol and btol. 2 damp is zero. The system A*x = b is probably\n not compatible. A least-squares solution has\n been obtained that is sufficiently accurate,\n given the value of atol. 3 damp is nonzero. A damped least-squares\n solution has been obtained that is sufficiently\n accurate, given the value of atol. 4 An estimate of cond(Abar) has exceeded conlim.\n The system A*x = b appears to be ill-conditioned,\n or there could be an error in Aprod1 or Aprod2. 5 The iteration limit itnlim was reached. integer(kind=ip), intent(out) :: itn The number of iterations performed. real(kind=dp), intent(out) :: normA An estimate of the Frobenius norm of Abar.\nThis is the square-root of the sum of squares\nof the elements of Abar.\nIf damp is small and the columns of A\nhave all been scaled to have length 1.0,\nnormA should increase to roughly sqrt(n).\nA radically different value for normA may\nindicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: condA An estimate of cond(Abar), the condition\nnumber of Abar. A very high value of condA\nmay again indicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: normr An estimate of the final value of norm(rbar),\nthe function being minimized (see notation\nabove). This will be small if A*x = b has\na solution. real(kind=dp), intent(out) :: normAr An estimate of the final value of norm( Abar'*rbar ) , the norm of\nthe residual for the normal equations.\nThis should be small in all cases. (normAr\nwill often be smaller than the true value\ncomputed from the output vector x.) real(kind=dp), intent(out) :: normx An estimate of norm(x) for the final solution x. Called by proc~~lsmr~~CalledByGraph proc~lsmr lsmrModule::lsmr proc~lsmr_ez lsmrModule::lsmr_ez proc~lsmr_ez->proc~lsmr Help Graph Key Nodes of different colours represent the following: Graph Key Subroutine Subroutine Function Function Interface Interface Type Bound Procedure Type Bound Procedure Unknown Procedure Type Unknown Procedure Type Program Program This Page's Entity This Page's Entity Solid arrows point from a procedure to one which it calls. Dashed \narrows point from an interface to procedures which implement that interface.\nThis could include the module procedures in a generic interface or the\nimplementation in a submodule of an interface in a parent module. Source Code subroutine lsmr ( m , n , Aprod1 , Aprod2 , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer ( ip ), intent ( in ) :: m !! the number of rows in A. integer ( ip ), intent ( in ) :: n !! the number of columns in A. integer ( ip ), intent ( in ) :: itnlim !! An upper limit on the number of iterations. !! Suggested value: !! !! * `itnlim = n/2` for well-conditioned systems !! with clustered singular values, !! * `itnlim = 4*n` otherwise. integer ( ip ), intent ( in ) :: localSize !! No. of vectors for local reorthogonalization: !! !! * 0 No reorthogonalization is performed. !! * >0 This many n-vectors \"v\" (the most recent ones) !! are saved for reorthogonalizing the next v. !! !! localSize need not be more than min(m,n). !! At most min(m,n) vectors will be allocated. integer ( ip ), intent ( in ) :: nout !! File number for printed output. If positive, !! a summary will be printed on file nout. integer ( ip ), intent ( out ) :: istop !! An integer giving the reason for termination: !! !! * 0 x = 0 is the exact solution. !! No iterations were performed. !! * 1 The equations `A*x = b` are probably compatible. !! `Norm(A*x - b)` is sufficiently small, given the !! values of atol and btol. !! * 2 damp is zero. The system `A*x = b` is probably !! not compatible. A least-squares solution has !! been obtained that is sufficiently accurate, !! given the value of atol. !! * 3 damp is nonzero. A damped least-squares !! solution has been obtained that is sufficiently !! accurate, given the value of atol. !! * 4 An estimate of cond(Abar) has exceeded conlim. !! The system `A*x = b` appears to be ill-conditioned, !! or there could be an error in Aprod1 or Aprod2. !! * 5 The iteration limit itnlim was reached. integer ( ip ), intent ( out ) :: itn !! The number of iterations performed. real ( dp ), intent ( in ) :: b ( m ) !! The rhs vector `b`. real ( dp ), intent ( out ) :: x ( n ) !! Returns the computed solution `x`. real ( dp ), intent ( in ) :: atol !! An estimate of the relative error in the data !! defining the matrix A. For example, if A is !! accurate to about 6 digits, set atol = 1.0e-6. real ( dp ), intent ( in ) :: btol !! An estimate of the relative error in the data !! defining the rhs b. For example, if b is !! accurate to about 6 digits, set btol = 1.0e-6. real ( dp ), intent ( in ) :: conlim !! An upper limit on cond(Abar), the apparent !! condition number of the matrix Abar. !! Iterations will be terminated if a computed !! estimate of cond(Abar) exceeds conlim. !! This is intended to prevent certain small or !! zero singular values of A or Abar from !! coming into effect and causing unwanted growth !! in the computed solution. !! !! conlim and damp may be used separately or !! together to regularize ill-conditioned systems. !! !! Normally, conlim should be in the range !! 1000 to 1/eps. !! Suggested value: !! !! * `conlim = 1/(100*eps)` for compatible systems, !! * `conlim = 1/(10*sqrt(eps))` for least squares. real ( dp ), intent ( in ) :: damp !! The damping parameter for problem 3 above. !! (damp should be 0.0 for problems 1 and 2.) !! If the system `A*x = b` is incompatible, values !! of damp in the range 0 to `sqrt(eps)*norm(A)` !! will probably have a negligible effect. !! Larger values of damp will tend to decrease !! the norm of x and reduce the number of !! iterations required by LSMR. !! !! The work per iteration and the storage needed !! by LSMR are the same for all values of damp. real ( dp ), intent ( out ) :: normA !! An estimate of the Frobenius norm of Abar. !! This is the square-root of the sum of squares !! of the elements of Abar. !! If damp is small and the columns of A !! have all been scaled to have length 1.0, !! normA should increase to roughly sqrt(n). !! A radically different value for normA may !! indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: condA !! An estimate of cond(Abar), the condition !! number of Abar. A very high value of condA !! may again indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: normr !! An estimate of the final value of norm(rbar), !! the function being minimized (see notation !! above). This will be small if A*x = b has !! a solution. real ( dp ), intent ( out ) :: normAr !! An estimate of the final value of !! `norm( Abar'*rbar )`, the norm of !! the residual for the normal equations. !! This should be small in all cases. (normAr !! will often be smaller than the true value !! computed from the output vector x.) real ( dp ), intent ( out ) :: normx !! An estimate of norm(x) for the final solution x. procedure ( Aprod1_f ) :: Aprod1 !! See above. procedure ( Aprod2_f ) :: Aprod2 !! See above. ! Local arrays and variables real ( dp ) :: h ( n ), hbar ( n ), u ( m ), v ( n ), w ( n ), localV ( n , min ( localSize , m , n )) logical :: damped , localOrtho , localVQueueFull , prnt , show integer ( ip ) :: i , localOrthoCount , localOrthoLimit , localPointer , localVecs , & pcount , pfreq real ( dp ) :: alpha , alphabar , alphahat , & beta , betaacute , betacheck , betad , betadd , betahat , & normb , c , cbar , chat , ctildeold , ctol , & d , maxrbar , minrbar , normA2 , & rho , rhobar , rhobarold , rhodold , rhoold , rhotemp , & rhotildeold , rtol , s , sbar , shat , stildeold , & t1 , taud , tautildeold , test1 , test2 , test3 , & thetabar , thetanew , thetatilde , thetatildeold , & zeta , zetabar , zetaold ! Local constants real ( dp ), parameter :: zero = 0.0_dp , one = 1.0_dp character ( len =* ), parameter :: enter = ' Enter LSMR. ' character ( len =* ), parameter :: exitt = ' Exit LSMR. ' character ( len =* ), parameter :: msg ( 0 : 7 ) = & ( / 'The exact solution is x = 0 ' , & 'Ax - b is small enough, given atol, btol ' , & 'The least-squares solution is good enough, given atol' , & 'The estimate of cond(Abar) has exceeded conlim ' , & 'Ax - b is small enough for this machine ' , & 'The LS solution is good enough for this machine ' , & 'Cond(Abar) seems to be too large for this machine ' , & 'The iteration limit has been reached ' / ) !------------------------------------------------------------------- ! Initialize. localVecs = min ( localSize , m , n ) show = nout > 0 if ( show ) then write ( nout , 1000 ) enter , m , n , damp , atol , conlim , btol , itnlim , localVecs end if pfreq = 20 ! print frequency (for repeating the heading) pcount = 0 ! print counter damped = damp > zero ! !------------------------------------------------------------------- ! Set up the first vectors u and v for the bidiagonalization. ! These satisfy beta*u = b, alpha*v = A(transpose)*u. !------------------------------------------------------------------- u ( 1 : m ) = b ( 1 : m ) v ( 1 : n ) = zero x ( 1 : n ) = zero alpha = zero beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) call Aprod2 ( m , n , v , u ) ! v = A'*u alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) end if if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) w = v end if normAr = alpha * beta if ( normAr == zero ) go to 800 ! Initialization for local reorthogonalization. localOrtho = . false . if ( localVecs > 0 ) then localPointer = 1 localOrtho = . true . localVQueueFull = . false . localV (:, 1 ) = v end if ! Initialize variables for 1st iteration. itn = 0 zetabar = alpha * beta alphabar = alpha rho = 1 rhobar = 1 cbar = 1 sbar = 0 h = v hbar ( 1 : n ) = zero x ( 1 : n ) = zero ! Initialize variables for estimation of ||r||. betadd = beta betad = 0 rhodold = 1 tautildeold = 0 thetatilde = 0 zeta = 0 d = 0 ! Initialize variables for estimation of ||A|| and cond(A). normA2 = alpha ** 2 maxrbar = 0_dp minrbar = 1e+100_dp ! Items for use in stopping rules. normb = beta istop = 0 ctol = zero if ( conlim > zero ) ctol = one / conlim normr = beta ! Exit if b=0 or A'b = 0. normAr = alpha * beta if ( normAr == 0 ) then if ( show ) then write ( nout , '(a)' ) msg ( 1 ) end if return end if ! Heading for iteration log. if ( show ) then if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if test1 = one test2 = alpha / beta write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 end if !=================================================================== ! Main iteration loop. !=================================================================== do itn = itn + 1 !---------------------------------------------------------------- ! Perform the next step of the bidiagonalization to obtain the ! next beta, u, alpha, v. These satisfy ! beta*u = A*v - alpha*u, ! alpha*v = A'*u - beta*v. !---------------------------------------------------------------- u = ( - alpha ) * u ! call dscal (m,(- alpha), u, 1) call Aprod1 ( m , n , v , u ) ! u = u + A*v beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) if ( localOrtho ) then ! Store v into the circular buffer localV. call localVEnqueue ! Store old v for local reorthog'n of new v. end if v = ( - beta ) * v ! call dscal (n, (- beta), v, 1) call Aprod2 ( m , n , v , u ) ! v = v + A'*u if ( localOrtho ) then ! Perform local reorthogonalization of V. call localVOrtho ! Local-reorthogonalization of new v. end if alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) end if end if ! At this point, beta = beta_{k+1}, alpha = alpha_{k+1}. !---------------------------------------------------------------- ! Construct rotation Qhat_{k,2k+1}. alphahat = d2norm ( alphabar , damp ) chat = alphabar / alphahat shat = damp / alphahat ! Use a plane rotation (Q_i) to turn B_i to R_i. rhoold = rho rho = d2norm ( alphahat , beta ) c = alphahat / rho s = beta / rho thetanew = s * alpha alphabar = c * alpha ! Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar. rhobarold = rhobar zetaold = zeta thetabar = sbar * rho rhotemp = cbar * rho rhobar = d2norm ( cbar * rho , thetanew ) cbar = cbar * rho / rhobar sbar = thetanew / rhobar zeta = cbar * zetabar zetabar = - sbar * zetabar ! Update h, h_hat, x. hbar = h - ( thetabar * rho / ( rhoold * rhobarold )) * hbar x = x + ( zeta / ( rho * rhobar )) * hbar h = v - ( thetanew / rho ) * h ! Estimate ||r||. ! Apply rotation Qhat_{k,2k+1}. betaacute = chat * betadd betacheck = - shat * betadd ! Apply rotation Q_{k,k+1}. betahat = c * betaacute betadd = - s * betaacute ! Apply rotation Qtilde_{k-1}. ! betad = betad_{k-1} here. thetatildeold = thetatilde rhotildeold = d2norm ( rhodold , thetabar ) ctildeold = rhodold / rhotildeold stildeold = thetabar / rhotildeold thetatilde = stildeold * rhobar rhodold = ctildeold * rhobar betad = - stildeold * betad + ctildeold * betahat ! betad = betad_k here. ! rhodold = rhod_k here. tautildeold = ( zetaold - thetatildeold * tautildeold ) / rhotildeold taud = ( zeta - thetatilde * tautildeold ) / rhodold d = d + betacheck ** 2 normr = sqrt ( d + ( betad - taud ) ** 2 + betadd ** 2 ) ! Estimate ||A||. normA2 = normA2 + beta ** 2 normA = sqrt ( normA2 ) normA2 = normA2 + alpha ** 2 ! Estimate cond(A). maxrbar = max ( maxrbar , rhobarold ) if ( itn > 1 ) then minrbar = min ( minrbar , rhobarold ) end if condA = max ( maxrbar , rhotemp ) / min ( minrbar , rhotemp ) !---------------------------------------------------------------- ! Test for convergence. !---------------------------------------------------------------- ! Compute norms for convergence testing. normAr = abs ( zetabar ) normx = sqrt ( dot_product ( x , x ) ) ! dnrm2(n, x, 1) ! Now use these norms to estimate certain other quantities, ! some of which will be small near a solution. test1 = normr / normb test2 = normAr / ( normA * normr ) test3 = one / condA t1 = test1 / ( one + normA * normx / normb ) rtol = btol + atol * normA * normx / normb ! The following tests guard against extremely small values of ! atol, btol or ctol. (The user may have set any or all of ! the parameters atol, btol, conlim to 0.) ! The effect is equivalent to the normAl tests using ! atol = eps, btol = eps, conlim = 1/eps. if ( itn >= itnlim ) istop = 7 if ( one + test3 <= one ) istop = 6 if ( one + test2 <= one ) istop = 5 if ( one + t1 <= one ) istop = 4 ! Allow for tolerances set by the user. if ( test3 <= ctol ) istop = 3 if ( test2 <= atol ) istop = 2 if ( test1 <= rtol ) istop = 1 !---------------------------------------------------------------- ! See if it is time to print something. !---------------------------------------------------------------- prnt = . false . if ( show ) then if ( n <= 40 ) prnt = . true . if ( itn <= 10 ) prnt = . true . if ( itn >= itnlim - 10 ) prnt = . true . if ( mod ( itn , 10 ) == 0 ) prnt = . true . if ( test3 <= 1.1 * ctol ) prnt = . true . if ( test2 <= 1.1 * atol ) prnt = . true . if ( test1 <= 1.1 * rtol ) prnt = . true . if ( istop /= 0 ) prnt = . true . if ( prnt ) then ! Print a line for this iteration if ( pcount >= pfreq ) then ! Print a heading first pcount = 0 if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if end if pcount = pcount + 1 write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 , normA , condA end if end if if ( istop /= 0 ) exit end do !=================================================================== ! End of iteration loop. !=================================================================== ! Come here if normAr = 0, or if normal exit. 800 if ( show ) then ! Print the stopping condition. write ( nout , 2000 ) & exitt , istop , itn , & exitt , normA , condA , & exitt , normb , normx , & exitt , normr , normAr write ( nout , 3000 ) & exitt , msg ( istop ) end if return 1000 format ( // a , ' Least-squares solution of Ax = b' & / ' The matrix A has' , i7 , ' rows and' , i7 , ' columns' & / ' damp =' , es22 . 14 & / ' atol =' , es10 . 2 , 15 x , 'conlim =' , es10 . 2 & / ' btol =' , es10 . 2 , 15 x , 'itnlim =' , i10 & / ' localSize (no. of vectors for local reorthogonalization) =' , i7 ) 1200 format ( / \" Itn x(1) norm r A'r \" , & ' Compatible LS norm A cond A' ) 1300 format ( / \" Itn x(1) norm rbar Abar'rbar\" , & ' Compatible LS norm Abar cond Abar' ) 1500 format ( i6 , 2 es17 . 9 , 5 es10 . 2 ) 2000 format ( / a , 5 x , 'istop =' , i2 , 15 x , 'itn =' , i8 & / a , 5 x , 'normA =' , es12 . 5 , 5 x , 'condA =' , es12 . 5 & / a , 5 x , 'normb =' , es12 . 5 , 5 x , 'normx =' , es12 . 5 & / a , 5 x , 'normr =' , es12 . 5 , 5 x , 'normAr =' , es12 . 5 ) 3000 format ( a , 5 x , a ) contains function d2norm ( a , b ) real ( dp ) :: d2norm real ( dp ), intent ( in ) :: a , b !------------------------------------------------------------------- ! d2norm returns sqrt( a**2 + b**2 ) ! with precautions to avoid overflow. ! ! 21 Mar 1990: First version. ! 17 Sep 2007: Fortran 90 version. ! 24 Oct 2007: User real(dp) instead of compiler option -r8. !------------------------------------------------------------------- intrinsic :: abs , sqrt real ( dp ) :: scale real ( dp ), parameter :: zero = 0.0_dp scale = abs ( a ) + abs ( b ) if ( scale == zero ) then d2norm = zero else d2norm = scale * sqrt (( a / scale ) ** 2 + ( b / scale ) ** 2 ) end if end function d2norm !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVEnqueue ! Store v into the circular buffer localV. if ( localPointer < localVecs ) then localPointer = localPointer + 1 else localPointer = 1 localVQueueFull = . true . end if localV (:, localPointer ) = v end subroutine localVEnqueue !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVOrtho ! Perform local reorthogonalization of current v. real ( dp ) :: d if ( localVQueueFull ) then localOrthoLimit = localVecs else localOrthoLimit = localPointer end if do localOrthoCount = 1 , localOrthoLimit d = dot_product ( v , localV (:, localOrthoCount )) v = v - d * localV (:, localOrthoCount ) end do end subroutine localVOrtho end subroutine lsmr","tags":"","loc":"proc/lsmr.html"},{"title":"lsmr_ez – LSMR","text":"public subroutine lsmr_ez(m, n, irow, icol, a, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) Easy interface to lsmr . Instead of specifying the Aprod1 , Aprod2 functions,\n the sparsity pattern ( irow , icol ) and nonzero elemenets\n of a are input. History JW : 1/24/2024 : created. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n integer, intent(in), dimension(:) :: irow row indices of nonzero elements of A integer, intent(in), dimension(:) :: icol column indices of nonzero elements of A real(kind=dp), intent(in), dimension(:) :: a nonzero elements of A real(kind=dp), intent(in) :: b (m) real(kind=dp), intent(in) :: damp real(kind=dp), intent(in) :: atol real(kind=dp), intent(in) :: btol real(kind=dp), intent(in) :: conlim integer(kind=ip), intent(in) :: itnlim integer(kind=ip), intent(in) :: localSize integer(kind=ip), intent(in) :: nout real(kind=dp), intent(out) :: x (n) integer(kind=ip), intent(out) :: istop integer(kind=ip), intent(out) :: itn real(kind=dp), intent(out) :: normA real(kind=dp), intent(out) :: condA real(kind=dp), intent(out) :: normr real(kind=dp), intent(out) :: normAr real(kind=dp), intent(out) :: normx Calls proc~~lsmr_ez~~CallsGraph proc~lsmr_ez lsmrModule::lsmr_ez proc~lsmr lsmrModule::lsmr proc~lsmr_ez->proc~lsmr Help Graph Key Nodes of different colours represent the following: Graph Key Subroutine Subroutine Function Function Interface Interface Type Bound Procedure Type Bound Procedure Unknown Procedure Type Unknown Procedure Type Program Program This Page's Entity This Page's Entity Solid arrows point from a procedure to one which it calls. Dashed \narrows point from an interface to procedures which implement that interface.\nThis could include the module procedures in a generic interface or the\nimplementation in a submodule of an interface in a parent module. Source Code subroutine lsmr_ez ( m , n , irow , icol , a , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer , dimension (:), intent ( in ) :: irow !! row indices of nonzero elements of `A` integer , dimension (:), intent ( in ) :: icol !! column indices of nonzero elements of `A` real ( dp ), dimension (:), intent ( in ) :: a !! nonzero elements of `A` integer ( ip ), intent ( in ) :: m , n , itnlim , localSize , nout integer ( ip ), intent ( out ) :: istop , itn real ( dp ), intent ( in ) :: b ( m ) real ( dp ), intent ( out ) :: x ( n ) real ( dp ), intent ( in ) :: atol , btol , conlim , damp real ( dp ), intent ( out ) :: normA , condA , normr , normAr , normx integer ( ip ) :: num_nonzero_elements !! number of nonzero elements in the matrix real ( dp ), dimension (:), allocatable :: Ax , Aty !! temp arrays if ( size ( irow ) == size ( icol ) . and . size ( irow ) == size ( a )) then num_nonzero_elements = size ( irow ) allocate ( Ax ( m )) allocate ( Aty ( n )) call lsmr ( m , n , Aprod1_ez , Aprod2_ez , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) else error stop 'inconsistent sizes of input arrays irow, icol, a' end if contains ! see code from LSQR subroutine Aprod1_ez ( m , n , x , y ) !! y := y + A*x integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index ! A*x: Ax = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Ax ( r ) = Ax ( r ) + a ( i ) * x ( c ) end do y = y + Ax end subroutine Aprod1_ez subroutine Aprod2_ez ( m , n , x , y ) !! x := x + A'*y integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index Aty = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Aty ( c ) = Aty ( c ) + a ( i ) * y ( r ) end do x = x + Aty end subroutine Aprod2_ez end subroutine lsmr_ez","tags":"","loc":"proc/lsmr_ez.html"},{"title":"ddot – LSMR","text":"public function ddot(n, dx, incx, dy, incy) Dot product of two vectors. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vectors. real(kind=wp), intent(in) :: dx (*) the first vector integer(kind=ip), intent(in) :: incx the increment between successive entries in DX. real(kind=wp), intent(in) :: dy (*) the second vector integer(kind=ip), intent(in) :: incy the increment between successive entries in DY. Return Value real(kind=wp) Source Code function ddot ( n , dx , incx , dy , incy ) implicit none integer , intent ( in ) :: n , incx , incy double precision , intent ( in ) :: dx ( * ), dy ( * ) double precision :: ddot end function ddot","tags":"","loc":"proc/ddot.html"},{"title":"dnrm2 – LSMR","text":"public function dnrm2(n, x, incx) The euclidean norm of a vector sqrt ( X' * X ) . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n real(kind=wp), intent(in) :: x (*) integer(kind=ip), intent(in) :: incx Return Value real(kind=wp) Source Code function dnrm2 ( n , dx , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: dx ( * ) double precision :: dnrm2 end function dnrm2","tags":"","loc":"proc/dnrm2.html"},{"title":"dcopy – LSMR","text":"public subroutine dcopy(n, dx, incx, dy, incy) Copies a vector X to a vector Y. Arguments Type Intent Optional Attributes Name integer(kind=ip) :: n real(kind=wp) :: dx (*) integer(kind=ip) :: incx real(kind=wp) :: dy (*) integer(kind=ip) :: incy Source Code subroutine dcopy ( n , dx , incx , dy , incy ) implicit none double precision dx ( * ), dy ( * ) integer i , incx , incy , ix , iy , m , n end subroutine dcopy","tags":"","loc":"proc/dcopy.html"},{"title":"dscal – LSMR","text":"public subroutine dscal(n, sa, x, incx) Scales a vector by a constant. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vector. real(kind=wp), intent(in) :: sa the multiplier. real(kind=wp), intent(inout) :: x (*) the vector to be scaled. integer(kind=ip), intent(in) :: incx the increment between successive entries of X. Source Code subroutine dscal ( n , sa , x , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: sa double precision , intent ( inout ) :: x ( * ) end subroutine dscal","tags":"","loc":"proc/dscal.html"},{"title":"lsmrModule – LSMR","text":"LSMR solves Ax = b or min ||Ax - b|| with or without damping,\n using the iterative algorithm of David Fong and Michael Saunders. Authors David Fong & Michael Saunders, Systems Optimization Laboratory (SOL) See also http://www.stanford.edu/group/SOL/software/lsmr.html History 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m. 07 Sep 2010: Local reorthogonalization now works (localSize > 0). 28 Jan 2014: In lsmrDataModule.f90:\n ip added for integer(ip) declarations.\n dnrm2 and dscal coded directly\n (no longer use lsmrblasInterface.f90 or lsmrblas.f90). Uses lsmrDataModule module~~lsmrmodule~~UsesGraph module~lsmrmodule lsmrModule module~lsmrdatamodule lsmrDataModule module~lsmrmodule->module~lsmrdatamodule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Abstract Interfaces abstract interface private subroutine Aprod1_f(m, n, x, y) y := y + A*x Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(in) :: x (n) real(kind=dp), intent(inout) :: y (m) abstract interface private subroutine Aprod2_f(m, n, x, y) x := x + A'*y Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n real(kind=dp), intent(inout) :: x (n) real(kind=dp), intent(in) :: y (m) Subroutines public subroutine lsmr (m, n, Aprod1, Aprod2, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) LSMR finds a solution x to the following problems: Read more… Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m the number of rows in A. integer(kind=ip), intent(in) :: n the number of columns in A. procedure( Aprod1_f ) :: Aprod1 See above. procedure( Aprod2_f ) :: Aprod2 See above. real(kind=dp), intent(in) :: b (m) The rhs vector b . real(kind=dp), intent(in) :: damp The damping parameter for problem 3 above.\n(damp should be 0.0 for problems 1 and 2.)\nIf the system A*x = b is incompatible, values\nof damp in the range 0 to sqrt(eps)*norm(A) will probably have a negligible effect.\nLarger values of damp will tend to decrease\nthe norm of x and reduce the number of\niterations required by LSMR. Read more… real(kind=dp), intent(in) :: atol An estimate of the relative error in the data\ndefining the matrix A. For example, if A is\naccurate to about 6 digits, set atol = 1.0e-6. real(kind=dp), intent(in) :: btol An estimate of the relative error in the data\ndefining the rhs b. For example, if b is\naccurate to about 6 digits, set btol = 1.0e-6. real(kind=dp), intent(in) :: conlim An upper limit on cond(Abar), the apparent\ncondition number of the matrix Abar.\nIterations will be terminated if a computed\nestimate of cond(Abar) exceeds conlim.\nThis is intended to prevent certain small or\nzero singular values of A or Abar from\ncoming into effect and causing unwanted growth\nin the computed solution. Read more… integer(kind=ip), intent(in) :: itnlim An upper limit on the number of iterations.\nSuggested value: Read more… integer(kind=ip), intent(in) :: localSize No. of vectors for local reorthogonalization: Read more… integer(kind=ip), intent(in) :: nout File number for printed output. If positive,\na summary will be printed on file nout. real(kind=dp), intent(out) :: x (n) Returns the computed solution x . integer(kind=ip), intent(out) :: istop An integer giving the reason for termination: Read more… integer(kind=ip), intent(out) :: itn The number of iterations performed. real(kind=dp), intent(out) :: normA An estimate of the Frobenius norm of Abar.\nThis is the square-root of the sum of squares\nof the elements of Abar.\nIf damp is small and the columns of A\nhave all been scaled to have length 1.0,\nnormA should increase to roughly sqrt(n).\nA radically different value for normA may\nindicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: condA An estimate of cond(Abar), the condition\nnumber of Abar. A very high value of condA\nmay again indicate an error in Aprod1 or Aprod2. real(kind=dp), intent(out) :: normr An estimate of the final value of norm(rbar),\nthe function being minimized (see notation\nabove). This will be small if A*x = b has\na solution. real(kind=dp), intent(out) :: normAr An estimate of the final value of norm( Abar'*rbar ) , the norm of\nthe residual for the normal equations.\nThis should be small in all cases. (normAr\nwill often be smaller than the true value\ncomputed from the output vector x.) real(kind=dp), intent(out) :: normx An estimate of norm(x) for the final solution x. public subroutine lsmr_ez (m, n, irow, icol, a, b, damp, atol, btol, conlim, itnlim, localSize, nout, x, istop, itn, normA, condA, normr, normAr, normx) Easy interface to lsmr . Read more… Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: m integer(kind=ip), intent(in) :: n integer, intent(in), dimension(:) :: irow row indices of nonzero elements of A integer, intent(in), dimension(:) :: icol column indices of nonzero elements of A real(kind=dp), intent(in), dimension(:) :: a nonzero elements of A real(kind=dp), intent(in) :: b (m) real(kind=dp), intent(in) :: damp real(kind=dp), intent(in) :: atol real(kind=dp), intent(in) :: btol real(kind=dp), intent(in) :: conlim integer(kind=ip), intent(in) :: itnlim integer(kind=ip), intent(in) :: localSize integer(kind=ip), intent(in) :: nout real(kind=dp), intent(out) :: x (n) integer(kind=ip), intent(out) :: istop integer(kind=ip), intent(out) :: itn real(kind=dp), intent(out) :: normA real(kind=dp), intent(out) :: condA real(kind=dp), intent(out) :: normr real(kind=dp), intent(out) :: normAr real(kind=dp), intent(out) :: normx","tags":"","loc":"module/lsmrmodule.html"},{"title":"lsmrblas – LSMR","text":"This file contains the following BLAS routines dcopy , ddot , dnrm2 , dscal required by subroutines lsmr and acheck . References Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart,\n LINPACK User's Guide,\n SIAM, 1979,\n ISBN13: 978-0-898711-72-1,\n LC: QA214.L56. Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh,\n Algorithm 539,\n Basic Linear Algebra Subprograms for Fortran Usage,\n ACM Transactions on Mathematical Software,\n Volume 5, Number 3, September 1979, pages 308-323. Note This version has not been linked with an external BLAS library,\n and uses the BLAS functions in this module. Uses lsmrDataModule module~~lsmrblas~~UsesGraph module~lsmrblas lsmrblas module~lsmrdatamodule lsmrDataModule module~lsmrblas->module~lsmrdatamodule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Functions public function ddot (n, dx, incx, dy, incy) Dot product of two vectors. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vectors. real(kind=wp), intent(in) :: dx (*) the first vector integer(kind=ip), intent(in) :: incx the increment between successive entries in DX. real(kind=wp), intent(in) :: dy (*) the second vector integer(kind=ip), intent(in) :: incy the increment between successive entries in DY. Return Value real(kind=wp) public function dnrm2 (n, x, incx) The euclidean norm of a vector sqrt ( X' * X ) . Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n real(kind=wp), intent(in) :: x (*) integer(kind=ip), intent(in) :: incx Return Value real(kind=wp) Subroutines public subroutine dcopy (n, dx, incx, dy, incy) Copies a vector X to a vector Y. Arguments Type Intent Optional Attributes Name integer(kind=ip) :: n real(kind=wp) :: dx (*) integer(kind=ip) :: incx real(kind=wp) :: dy (*) integer(kind=ip) :: incy public subroutine dscal (n, sa, x, incx) Scales a vector by a constant. Arguments Type Intent Optional Attributes Name integer(kind=ip), intent(in) :: n the number of entries in the vector. real(kind=wp), intent(in) :: sa the multiplier. real(kind=wp), intent(inout) :: x (*) the vector to be scaled. integer(kind=ip), intent(in) :: incx the increment between successive entries of X.","tags":"","loc":"module/lsmrblas.html"},{"title":"lsmrDataModule – LSMR","text":"Defines integer(ip) and real(wp) kind parameters. Notes JW: 1/24/2024 : The original version of this file was rewritten. Uses iso_fortran_env module~~lsmrdatamodule~~UsesGraph module~lsmrdatamodule lsmrDataModule iso_fortran_env iso_fortran_env module~lsmrdatamodule->iso_fortran_env Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Used by module~~lsmrdatamodule~~UsedByGraph module~lsmrdatamodule lsmrDataModule module~lsmrblas lsmrblas module~lsmrblas->module~lsmrdatamodule module~lsmrmodule lsmrModule module~lsmrmodule->module~lsmrdatamodule Help Graph Key Nodes of different colours represent the following: Graph Key Module Module Submodule Submodule Subroutine Subroutine Function Function Program Program This Page's Entity This Page's Entity Solid arrows point from a submodule to the (sub)module which it is\ndescended from. Dashed arrows point from a module or program unit to \nmodules which it uses. Variables Type Visibility Attributes Name Initial integer, public, parameter :: lsmr_wp = real64 real kind used by this module [8 bytes] integer, public, parameter :: lsmr_ip = int32 integer kind used by this module [4 bytes]","tags":"","loc":"module/lsmrdatamodule.html"},{"title":"lsmrModule.f90 – LSMR","text":"This file depends on sourcefile~~lsmrmodule.f90~~EfferentGraph sourcefile~lsmrmodule.f90 lsmrModule.f90 sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrmodule.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! LSMR solves `Ax = b` or `min ||Ax - b||` with or without damping, ! using the iterative algorithm of David Fong and Michael Saunders. ! !### Authors ! * David Fong & Michael Saunders, Systems Optimization Laboratory (SOL) ! !### See also ! * http://www.stanford.edu/group/SOL/software/lsmr.html ! !### History ! * 17 Jul 2010: F90 LSMR derived from F90 LSQR and lsqr.m. ! * 07 Sep 2010: Local reorthogonalization now works (localSize > 0). ! * 28 Jan 2014: In lsmrDataModule.f90: ! ip added for integer(ip) declarations. ! dnrm2 and dscal coded directly ! (no longer use lsmrblasInterface.f90 or lsmrblas.f90). module lsmrModule use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none private public :: lsmr public :: lsmr_ez abstract interface subroutine Aprod1_f ( m , n , x , y ) !! y := y + A*x use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) end subroutine Aprod1_f subroutine Aprod2_f ( m , n , x , y ) !! x := x + A'*y use lsmrDataModule , only : ip => lsmr_ip , dp => lsmr_wp implicit none integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) end subroutine Aprod2_f end interface contains !***************************************************************************** !> ! LSMR finds a solution x to the following problems: ! ! 1. Unsymmetric equations: Solve `A*x = b` ! 2. Linear least squares: Solve `A*x = b` in the least-squares sense ! 3. Damped least squares: Solve !``` ! ( A )*x = ( b ) ! ( damp*I ) ( 0 ) !``` ! in the least-squares sense ! ! where A is a matrix with m rows and n columns, b is an m-vector, ! and damp is a scalar. (All quantities are real.) ! The matrix A is treated as a linear operator. It is accessed ! by means of subroutine calls with the following purpose: ! ! * `call Aprod1(m,n,x,y)` must compute `y = y + A*x ` without altering `x`. ! * `call Aprod2(m,n,x,y)` must compute `x = x + A'*y` without altering `y`. ! ! LSMR uses an iterative method to approximate the solution. ! The number of iterations required to reach a certain accuracy ! depends strongly on the scaling of the problem. Poor scaling of ! the rows or columns of A should therefore be avoided where ! possible. ! ! For example, in problem 1 the solution is unaltered by ! row-scaling. If a row of A is very small or large compared to ! the other rows of A, the corresponding row of ( A b ) should be ! scaled up or down. ! ! In problems 1 and 2, the solution x is easily recovered ! following column-scaling. Unless better information is known, ! the nonzero columns of A should be scaled so that they all have ! the same Euclidean norm (e.g., 1.0). ! ! In problem 3, there is no freedom to re-scale if damp is ! nonzero. However, the value of damp should be assigned only ! after attention has been paid to the scaling of A. ! ! The parameter damp is intended to help regularize ! ill-conditioned systems, by preventing the true solution from ! being very large. Another aid to regularization is provided by ! the parameter condA, which may be used to terminate iterations ! before the computed solution becomes very large. ! ! Note that x is not an input parameter. ! If some initial estimate x0 is known and if damp = 0, ! one could proceed as follows: ! ! 1. Compute a residual vector `r0 = b - A*x0`. ! 2. Use LSMR to solve the system `A*dx = r0`. ! 3. Add the correction `dx` to obtain a final solution `x = x0 + dx`. ! ! This requires that x0 be available before and after the call ! to LSMR. To judge the benefits, suppose LSMR takes k1 iterations ! to solve `A*x = b` and k2 iterations to solve `A*dx = r0`. ! If x0 is \"good\", norm(r0) will be smaller than norm(b). ! If the same stopping tolerances atol and btol are used for each ! system, k1 and k2 will be similar, but the final solution `x0 + dx` ! should be more accurate. The only way to reduce the total work ! is to use a larger stopping tolerance for the second system. ! If some value btol is suitable for `A*x = b`, the larger value ! `btol*norm(b)/norm(r0)` should be suitable for `A*dx = r0`. ! ! Preconditioning is another way to reduce the number of iterations. ! If it is possible to solve a related system `M*x = b` efficiently, ! where M approximates A in some helpful way ! (e.g. M - A has low rank or its elements are small relative to ! those of A), LSMR may converge more rapidly on the system ! `A*M(inverse)*z = b`, ! after which `x` can be recovered by solving `M*x = z`. ! ! NOTE: If A is symmetric, LSMR should not be used! ! Alternatives are the symmetric conjugate-gradient method (CG) ! and/or SYMMLQ. ! SYMMLQ is an implementation of symmetric CG that applies to ! any symmetric A and will converge more rapidly than LSMR. ! If A is positive definite, there are other implementations of ! symmetric CG that require slightly less work per iteration ! than SYMMLQ (but will take the same number of iterations). ! ! Notation ! -------- ! The following quantities are used in discussing the subroutine ! parameters: !``` ! Abar = ( A ), bbar = (b) ! (damp*I) (0) ! ! r = b - A*x, rbar = bbar - Abar*x ! ! normr = sqrt( norm(r)**2 + damp**2 * norm(x)**2 ) ! = norm( rbar ) ! ! eps = the relative precision of floating-point arithmetic. ! On most machines, eps is about 1.0e-7 and 1.0e-16 ! in single and double precision respectively. ! We expect eps to be about 1e-16 always. !``` ! LSMR minimizes the function `normr` with respect to `x`. ! !### Precision ! ! The number of iterations required by LSMR will decrease ! if the computation is performed in higher precision. ! !### Reference ! * http://www.stanford.edu/group/SOL/software/lsmr.html ! !### LSMR development: ! * 21 Sep 2007: Fortran 90 version of LSQR implemented. ! Aprod1, Aprod2 implemented via f90 interface. ! * 17 Jul 2010: LSMR derived from LSQR and lsmr.m. ! * 07 Sep 2010: Local reorthogonalization now working. ! * 02 May 2014: With damp>0, istop=2 was incorrectly set to istop=3 ! (so incorrect stopping message was printed). Fixed. ! !@note Any or all of `atol`, `btol`, `conlim` may be set to zero. ! The effect will be the same as the values `eps`, `eps`, `1/eps`. subroutine lsmr ( m , n , Aprod1 , Aprod2 , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer ( ip ), intent ( in ) :: m !! the number of rows in A. integer ( ip ), intent ( in ) :: n !! the number of columns in A. integer ( ip ), intent ( in ) :: itnlim !! An upper limit on the number of iterations. !! Suggested value: !! !! * `itnlim = n/2` for well-conditioned systems !! with clustered singular values, !! * `itnlim = 4*n` otherwise. integer ( ip ), intent ( in ) :: localSize !! No. of vectors for local reorthogonalization: !! !! * 0 No reorthogonalization is performed. !! * >0 This many n-vectors \"v\" (the most recent ones) !! are saved for reorthogonalizing the next v. !! !! localSize need not be more than min(m,n). !! At most min(m,n) vectors will be allocated. integer ( ip ), intent ( in ) :: nout !! File number for printed output. If positive, !! a summary will be printed on file nout. integer ( ip ), intent ( out ) :: istop !! An integer giving the reason for termination: !! !! * 0 x = 0 is the exact solution. !! No iterations were performed. !! * 1 The equations `A*x = b` are probably compatible. !! `Norm(A*x - b)` is sufficiently small, given the !! values of atol and btol. !! * 2 damp is zero. The system `A*x = b` is probably !! not compatible. A least-squares solution has !! been obtained that is sufficiently accurate, !! given the value of atol. !! * 3 damp is nonzero. A damped least-squares !! solution has been obtained that is sufficiently !! accurate, given the value of atol. !! * 4 An estimate of cond(Abar) has exceeded conlim. !! The system `A*x = b` appears to be ill-conditioned, !! or there could be an error in Aprod1 or Aprod2. !! * 5 The iteration limit itnlim was reached. integer ( ip ), intent ( out ) :: itn !! The number of iterations performed. real ( dp ), intent ( in ) :: b ( m ) !! The rhs vector `b`. real ( dp ), intent ( out ) :: x ( n ) !! Returns the computed solution `x`. real ( dp ), intent ( in ) :: atol !! An estimate of the relative error in the data !! defining the matrix A. For example, if A is !! accurate to about 6 digits, set atol = 1.0e-6. real ( dp ), intent ( in ) :: btol !! An estimate of the relative error in the data !! defining the rhs b. For example, if b is !! accurate to about 6 digits, set btol = 1.0e-6. real ( dp ), intent ( in ) :: conlim !! An upper limit on cond(Abar), the apparent !! condition number of the matrix Abar. !! Iterations will be terminated if a computed !! estimate of cond(Abar) exceeds conlim. !! This is intended to prevent certain small or !! zero singular values of A or Abar from !! coming into effect and causing unwanted growth !! in the computed solution. !! !! conlim and damp may be used separately or !! together to regularize ill-conditioned systems. !! !! Normally, conlim should be in the range !! 1000 to 1/eps. !! Suggested value: !! !! * `conlim = 1/(100*eps)` for compatible systems, !! * `conlim = 1/(10*sqrt(eps))` for least squares. real ( dp ), intent ( in ) :: damp !! The damping parameter for problem 3 above. !! (damp should be 0.0 for problems 1 and 2.) !! If the system `A*x = b` is incompatible, values !! of damp in the range 0 to `sqrt(eps)*norm(A)` !! will probably have a negligible effect. !! Larger values of damp will tend to decrease !! the norm of x and reduce the number of !! iterations required by LSMR. !! !! The work per iteration and the storage needed !! by LSMR are the same for all values of damp. real ( dp ), intent ( out ) :: normA !! An estimate of the Frobenius norm of Abar. !! This is the square-root of the sum of squares !! of the elements of Abar. !! If damp is small and the columns of A !! have all been scaled to have length 1.0, !! normA should increase to roughly sqrt(n). !! A radically different value for normA may !! indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: condA !! An estimate of cond(Abar), the condition !! number of Abar. A very high value of condA !! may again indicate an error in Aprod1 or Aprod2. real ( dp ), intent ( out ) :: normr !! An estimate of the final value of norm(rbar), !! the function being minimized (see notation !! above). This will be small if A*x = b has !! a solution. real ( dp ), intent ( out ) :: normAr !! An estimate of the final value of !! `norm( Abar'*rbar )`, the norm of !! the residual for the normal equations. !! This should be small in all cases. (normAr !! will often be smaller than the true value !! computed from the output vector x.) real ( dp ), intent ( out ) :: normx !! An estimate of norm(x) for the final solution x. procedure ( Aprod1_f ) :: Aprod1 !! See above. procedure ( Aprod2_f ) :: Aprod2 !! See above. ! Local arrays and variables real ( dp ) :: h ( n ), hbar ( n ), u ( m ), v ( n ), w ( n ), localV ( n , min ( localSize , m , n )) logical :: damped , localOrtho , localVQueueFull , prnt , show integer ( ip ) :: i , localOrthoCount , localOrthoLimit , localPointer , localVecs , & pcount , pfreq real ( dp ) :: alpha , alphabar , alphahat , & beta , betaacute , betacheck , betad , betadd , betahat , & normb , c , cbar , chat , ctildeold , ctol , & d , maxrbar , minrbar , normA2 , & rho , rhobar , rhobarold , rhodold , rhoold , rhotemp , & rhotildeold , rtol , s , sbar , shat , stildeold , & t1 , taud , tautildeold , test1 , test2 , test3 , & thetabar , thetanew , thetatilde , thetatildeold , & zeta , zetabar , zetaold ! Local constants real ( dp ), parameter :: zero = 0.0_dp , one = 1.0_dp character ( len =* ), parameter :: enter = ' Enter LSMR. ' character ( len =* ), parameter :: exitt = ' Exit LSMR. ' character ( len =* ), parameter :: msg ( 0 : 7 ) = & ( / 'The exact solution is x = 0 ' , & 'Ax - b is small enough, given atol, btol ' , & 'The least-squares solution is good enough, given atol' , & 'The estimate of cond(Abar) has exceeded conlim ' , & 'Ax - b is small enough for this machine ' , & 'The LS solution is good enough for this machine ' , & 'Cond(Abar) seems to be too large for this machine ' , & 'The iteration limit has been reached ' / ) !------------------------------------------------------------------- ! Initialize. localVecs = min ( localSize , m , n ) show = nout > 0 if ( show ) then write ( nout , 1000 ) enter , m , n , damp , atol , conlim , btol , itnlim , localVecs end if pfreq = 20 ! print frequency (for repeating the heading) pcount = 0 ! print counter damped = damp > zero ! !------------------------------------------------------------------- ! Set up the first vectors u and v for the bidiagonalization. ! These satisfy beta*u = b, alpha*v = A(transpose)*u. !------------------------------------------------------------------- u ( 1 : m ) = b ( 1 : m ) v ( 1 : n ) = zero x ( 1 : n ) = zero alpha = zero beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) call Aprod2 ( m , n , v , u ) ! v = A'*u alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) end if if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) w = v end if normAr = alpha * beta if ( normAr == zero ) go to 800 ! Initialization for local reorthogonalization. localOrtho = . false . if ( localVecs > 0 ) then localPointer = 1 localOrtho = . true . localVQueueFull = . false . localV (:, 1 ) = v end if ! Initialize variables for 1st iteration. itn = 0 zetabar = alpha * beta alphabar = alpha rho = 1 rhobar = 1 cbar = 1 sbar = 0 h = v hbar ( 1 : n ) = zero x ( 1 : n ) = zero ! Initialize variables for estimation of ||r||. betadd = beta betad = 0 rhodold = 1 tautildeold = 0 thetatilde = 0 zeta = 0 d = 0 ! Initialize variables for estimation of ||A|| and cond(A). normA2 = alpha ** 2 maxrbar = 0_dp minrbar = 1e+100_dp ! Items for use in stopping rules. normb = beta istop = 0 ctol = zero if ( conlim > zero ) ctol = one / conlim normr = beta ! Exit if b=0 or A'b = 0. normAr = alpha * beta if ( normAr == 0 ) then if ( show ) then write ( nout , '(a)' ) msg ( 1 ) end if return end if ! Heading for iteration log. if ( show ) then if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if test1 = one test2 = alpha / beta write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 end if !=================================================================== ! Main iteration loop. !=================================================================== do itn = itn + 1 !---------------------------------------------------------------- ! Perform the next step of the bidiagonalization to obtain the ! next beta, u, alpha, v. These satisfy ! beta*u = A*v - alpha*u, ! alpha*v = A'*u - beta*v. !---------------------------------------------------------------- u = ( - alpha ) * u ! call dscal (m,(- alpha), u, 1) call Aprod1 ( m , n , v , u ) ! u = u + A*v beta = sqrt ( dot_product ( u , u ) ) ! dnrm2 (m, u, 1) if ( beta > zero ) then u = ( one / beta ) * u ! call dscal (m, (one/beta), u, 1) if ( localOrtho ) then ! Store v into the circular buffer localV. call localVEnqueue ! Store old v for local reorthog'n of new v. end if v = ( - beta ) * v ! call dscal (n, (- beta), v, 1) call Aprod2 ( m , n , v , u ) ! v = v + A'*u if ( localOrtho ) then ! Perform local reorthogonalization of V. call localVOrtho ! Local-reorthogonalization of new v. end if alpha = sqrt ( dot_product ( v , v ) ) ! dnrm2 (n, v, 1) if ( alpha > zero ) then v = ( one / alpha ) * v ! call dscal (n, (one/alpha), v, 1) end if end if ! At this point, beta = beta_{k+1}, alpha = alpha_{k+1}. !---------------------------------------------------------------- ! Construct rotation Qhat_{k,2k+1}. alphahat = d2norm ( alphabar , damp ) chat = alphabar / alphahat shat = damp / alphahat ! Use a plane rotation (Q_i) to turn B_i to R_i. rhoold = rho rho = d2norm ( alphahat , beta ) c = alphahat / rho s = beta / rho thetanew = s * alpha alphabar = c * alpha ! Use a plane rotation (Qbar_i) to turn R_i^T into R_i^bar. rhobarold = rhobar zetaold = zeta thetabar = sbar * rho rhotemp = cbar * rho rhobar = d2norm ( cbar * rho , thetanew ) cbar = cbar * rho / rhobar sbar = thetanew / rhobar zeta = cbar * zetabar zetabar = - sbar * zetabar ! Update h, h_hat, x. hbar = h - ( thetabar * rho / ( rhoold * rhobarold )) * hbar x = x + ( zeta / ( rho * rhobar )) * hbar h = v - ( thetanew / rho ) * h ! Estimate ||r||. ! Apply rotation Qhat_{k,2k+1}. betaacute = chat * betadd betacheck = - shat * betadd ! Apply rotation Q_{k,k+1}. betahat = c * betaacute betadd = - s * betaacute ! Apply rotation Qtilde_{k-1}. ! betad = betad_{k-1} here. thetatildeold = thetatilde rhotildeold = d2norm ( rhodold , thetabar ) ctildeold = rhodold / rhotildeold stildeold = thetabar / rhotildeold thetatilde = stildeold * rhobar rhodold = ctildeold * rhobar betad = - stildeold * betad + ctildeold * betahat ! betad = betad_k here. ! rhodold = rhod_k here. tautildeold = ( zetaold - thetatildeold * tautildeold ) / rhotildeold taud = ( zeta - thetatilde * tautildeold ) / rhodold d = d + betacheck ** 2 normr = sqrt ( d + ( betad - taud ) ** 2 + betadd ** 2 ) ! Estimate ||A||. normA2 = normA2 + beta ** 2 normA = sqrt ( normA2 ) normA2 = normA2 + alpha ** 2 ! Estimate cond(A). maxrbar = max ( maxrbar , rhobarold ) if ( itn > 1 ) then minrbar = min ( minrbar , rhobarold ) end if condA = max ( maxrbar , rhotemp ) / min ( minrbar , rhotemp ) !---------------------------------------------------------------- ! Test for convergence. !---------------------------------------------------------------- ! Compute norms for convergence testing. normAr = abs ( zetabar ) normx = sqrt ( dot_product ( x , x ) ) ! dnrm2(n, x, 1) ! Now use these norms to estimate certain other quantities, ! some of which will be small near a solution. test1 = normr / normb test2 = normAr / ( normA * normr ) test3 = one / condA t1 = test1 / ( one + normA * normx / normb ) rtol = btol + atol * normA * normx / normb ! The following tests guard against extremely small values of ! atol, btol or ctol. (The user may have set any or all of ! the parameters atol, btol, conlim to 0.) ! The effect is equivalent to the normAl tests using ! atol = eps, btol = eps, conlim = 1/eps. if ( itn >= itnlim ) istop = 7 if ( one + test3 <= one ) istop = 6 if ( one + test2 <= one ) istop = 5 if ( one + t1 <= one ) istop = 4 ! Allow for tolerances set by the user. if ( test3 <= ctol ) istop = 3 if ( test2 <= atol ) istop = 2 if ( test1 <= rtol ) istop = 1 !---------------------------------------------------------------- ! See if it is time to print something. !---------------------------------------------------------------- prnt = . false . if ( show ) then if ( n <= 40 ) prnt = . true . if ( itn <= 10 ) prnt = . true . if ( itn >= itnlim - 10 ) prnt = . true . if ( mod ( itn , 10 ) == 0 ) prnt = . true . if ( test3 <= 1.1 * ctol ) prnt = . true . if ( test2 <= 1.1 * atol ) prnt = . true . if ( test1 <= 1.1 * rtol ) prnt = . true . if ( istop /= 0 ) prnt = . true . if ( prnt ) then ! Print a line for this iteration if ( pcount >= pfreq ) then ! Print a heading first pcount = 0 if ( damped ) then write ( nout , 1300 ) else write ( nout , 1200 ) end if end if pcount = pcount + 1 write ( nout , 1500 ) itn , x ( 1 ), normr , normAr , test1 , test2 , normA , condA end if end if if ( istop /= 0 ) exit end do !=================================================================== ! End of iteration loop. !=================================================================== ! Come here if normAr = 0, or if normal exit. 800 if ( show ) then ! Print the stopping condition. write ( nout , 2000 ) & exitt , istop , itn , & exitt , normA , condA , & exitt , normb , normx , & exitt , normr , normAr write ( nout , 3000 ) & exitt , msg ( istop ) end if return 1000 format ( // a , ' Least-squares solution of Ax = b' & / ' The matrix A has' , i7 , ' rows and' , i7 , ' columns' & / ' damp =' , es22 . 14 & / ' atol =' , es10 . 2 , 15 x , 'conlim =' , es10 . 2 & / ' btol =' , es10 . 2 , 15 x , 'itnlim =' , i10 & / ' localSize (no. of vectors for local reorthogonalization) =' , i7 ) 1200 format ( / \" Itn x(1) norm r A'r \" , & ' Compatible LS norm A cond A' ) 1300 format ( / \" Itn x(1) norm rbar Abar'rbar\" , & ' Compatible LS norm Abar cond Abar' ) 1500 format ( i6 , 2 es17 . 9 , 5 es10 . 2 ) 2000 format ( / a , 5 x , 'istop =' , i2 , 15 x , 'itn =' , i8 & / a , 5 x , 'normA =' , es12 . 5 , 5 x , 'condA =' , es12 . 5 & / a , 5 x , 'normb =' , es12 . 5 , 5 x , 'normx =' , es12 . 5 & / a , 5 x , 'normr =' , es12 . 5 , 5 x , 'normAr =' , es12 . 5 ) 3000 format ( a , 5 x , a ) contains function d2norm ( a , b ) real ( dp ) :: d2norm real ( dp ), intent ( in ) :: a , b !------------------------------------------------------------------- ! d2norm returns sqrt( a**2 + b**2 ) ! with precautions to avoid overflow. ! ! 21 Mar 1990: First version. ! 17 Sep 2007: Fortran 90 version. ! 24 Oct 2007: User real(dp) instead of compiler option -r8. !------------------------------------------------------------------- intrinsic :: abs , sqrt real ( dp ) :: scale real ( dp ), parameter :: zero = 0.0_dp scale = abs ( a ) + abs ( b ) if ( scale == zero ) then d2norm = zero else d2norm = scale * sqrt (( a / scale ) ** 2 + ( b / scale ) ** 2 ) end if end function d2norm !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVEnqueue ! Store v into the circular buffer localV. if ( localPointer < localVecs ) then localPointer = localPointer + 1 else localPointer = 1 localVQueueFull = . true . end if localV (:, localPointer ) = v end subroutine localVEnqueue !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ subroutine localVOrtho ! Perform local reorthogonalization of current v. real ( dp ) :: d if ( localVQueueFull ) then localOrthoLimit = localVecs else localOrthoLimit = localPointer end if do localOrthoCount = 1 , localOrthoLimit d = dot_product ( v , localV (:, localOrthoCount )) v = v - d * localV (:, localOrthoCount ) end do end subroutine localVOrtho end subroutine lsmr !***************************************************************************** !> ! Easy interface to [[lsmr]]. ! ! Instead of specifying the `Aprod1`, `Aprod2` functions, ! the sparsity pattern (`irow`, `icol`) and nonzero elemenets ! of `a` are input. ! !### History ! * JW : 1/24/2024 : created. subroutine lsmr_ez ( m , n , irow , icol , a , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) integer , dimension (:), intent ( in ) :: irow !! row indices of nonzero elements of `A` integer , dimension (:), intent ( in ) :: icol !! column indices of nonzero elements of `A` real ( dp ), dimension (:), intent ( in ) :: a !! nonzero elements of `A` integer ( ip ), intent ( in ) :: m , n , itnlim , localSize , nout integer ( ip ), intent ( out ) :: istop , itn real ( dp ), intent ( in ) :: b ( m ) real ( dp ), intent ( out ) :: x ( n ) real ( dp ), intent ( in ) :: atol , btol , conlim , damp real ( dp ), intent ( out ) :: normA , condA , normr , normAr , normx integer ( ip ) :: num_nonzero_elements !! number of nonzero elements in the matrix real ( dp ), dimension (:), allocatable :: Ax , Aty !! temp arrays if ( size ( irow ) == size ( icol ) . and . size ( irow ) == size ( a )) then num_nonzero_elements = size ( irow ) allocate ( Ax ( m )) allocate ( Aty ( n )) call lsmr ( m , n , Aprod1_ez , Aprod2_ez , b , damp , & atol , btol , conlim , itnlim , localSize , nout , & x , istop , itn , normA , condA , normr , normAr , normx ) else error stop 'inconsistent sizes of input arrays irow, icol, a' end if contains ! see code from LSQR subroutine Aprod1_ez ( m , n , x , y ) !! y := y + A*x integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( in ) :: x ( n ) real ( dp ), intent ( inout ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index ! A*x: Ax = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Ax ( r ) = Ax ( r ) + a ( i ) * x ( c ) end do y = y + Ax end subroutine Aprod1_ez subroutine Aprod2_ez ( m , n , x , y ) !! x := x + A'*y integer ( ip ), intent ( in ) :: m , n real ( dp ), intent ( inout ) :: x ( n ) real ( dp ), intent ( in ) :: y ( m ) integer ( ip ) :: i !! counter integer ( ip ) :: r !! row index integer ( ip ) :: c !! column index Aty = 0.0_dp do i = 1 , num_nonzero_elements r = irow ( i ) c = icol ( i ) Aty ( c ) = Aty ( c ) + a ( i ) * y ( r ) end do x = x + Aty end subroutine Aprod2_ez end subroutine lsmr_ez !***************************************************************************** end module LSMRmodule","tags":"","loc":"sourcefile/lsmrmodule.f90.html"},{"title":"lsmrblas.F90 – LSMR","text":"This file depends on sourcefile~~lsmrblas.f90~~EfferentGraph sourcefile~lsmrblas.f90 lsmrblas.F90 sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrblas.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! This file contains the following BLAS routines ! [[dcopy]], [[ddot]], [[dnrm2]], [[dscal]] ! required by subroutines [[lsmr]] and [[acheck]]. ! !### References ! ! * Jack Dongarra, Jim Bunch, Cleve Moler, Pete Stewart, ! LINPACK User's Guide, ! SIAM, 1979, ! ISBN13: 978-0-898711-72-1, ! LC: QA214.L56. ! * Charles Lawson, Richard Hanson, David Kincaid, Fred Krogh, ! Algorithm 539, ! Basic Linear Algebra Subprograms for Fortran Usage, ! ACM Transactions on Mathematical Software, ! Volume 5, Number 3, September 1979, pages 308-323. ! #ifdef HAS_BLAS !@note This version has been linked with an external BLAS library, ! only the interfaces are defined in this module. #else !@note This version has not been linked with an external BLAS library, ! and uses the BLAS functions in this module. #endif module lsmrblas use lsmrDataModule , only : ip => lsmr_ip , wp => lsmr_wp implicit none private public :: ddot , dnrm2 , dscal , dcopy #ifdef HAS_BLAS ! get BLAS from an externally-linked library (double precision only) interface function ddot ( n , dx , incx , dy , incy ) implicit none integer , intent ( in ) :: n , incx , incy double precision , intent ( in ) :: dx ( * ), dy ( * ) double precision :: ddot end function ddot function dnrm2 ( n , dx , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: dx ( * ) double precision :: dnrm2 end function dnrm2 subroutine dscal ( n , sa , x , incx ) implicit none integer , intent ( in ) :: n , incx double precision , intent ( in ) :: sa double precision , intent ( inout ) :: x ( * ) end subroutine dscal subroutine dcopy ( n , dx , incx , dy , incy ) implicit none double precision dx ( * ), dy ( * ) integer i , incx , incy , ix , iy , m , n end subroutine dcopy end interface #else contains !***************************************************************************** !> ! Copies a vector X to a vector Y. subroutine dcopy ( n , dx , incx , dy , incy ) implicit none real ( wp ) dx ( * ), dy ( * ) integer ( ip ) i , incx , incy , ix , iy , m , n if ( n <= 0 ) then return end if if ( incx == 1 . and . incy == 1 ) then m = mod ( n , 7 ) if ( m /= 0 ) then dy ( 1 : m ) = dx ( 1 : m ) end if do i = m + 1 , n , 7 dy ( i ) = dx ( i ) dy ( i + 1 ) = dx ( i + 1 ) dy ( i + 2 ) = dx ( i + 2 ) dy ( i + 3 ) = dx ( i + 3 ) dy ( i + 4 ) = dx ( i + 4 ) dy ( i + 5 ) = dx ( i + 5 ) dy ( i + 6 ) = dx ( i + 6 ) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1 , n dy ( iy ) = dx ( ix ) ix = ix + incx iy = iy + incy end do end if end subroutine dcopy !***************************************************************************** !> ! Dot product of two vectors. real ( wp ) function ddot ( n , dx , incx , dy , incy ) integer ( ip ), intent ( in ) :: n !! the number of entries in the vectors. real ( wp ), intent ( in ) :: dx ( * ) !! the first vector integer ( ip ), intent ( in ) :: incx !! the increment between successive entries in DX. real ( wp ), intent ( in ) :: dy ( * ) !! the second vector integer ( ip ), intent ( in ) :: incy !! the increment between successive entries in DY. integer ( ip ) :: i , ix , iy , m real ( wp ) :: dtemp ddot = 0.0_wp dtemp = 0.0_wp if ( n <= 0 ) return if ( incx /= 1 . or . incy /= 1 ) then ! Code for unequal increments or equal increments ! not equal to 1. if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if if ( 0 <= incy ) then iy = 1 else iy = ( - n + 1 ) * incy + 1 end if do i = 1 , n dtemp = dtemp + dx ( ix ) * dy ( iy ) ix = ix + incx iy = iy + incy end do else ! Code for both increments equal to 1. m = mod ( n , 5 ) do i = 1 , m dtemp = dtemp + dx ( i ) * dy ( i ) end do do i = m + 1 , n , 5 dtemp = dtemp + dx ( i ) * dy ( i ) + dx ( i + 1 ) * dy ( i + 1 ) + dx ( i + 2 ) * dy ( i + 2 ) & + dx ( i + 3 ) * dy ( i + 3 ) + dx ( i + 4 ) * dy ( i + 4 ) end do end if ddot = dtemp end function ddot !***************************************************************************** !> ! The euclidean norm of a vector `sqrt ( X' * X )`. real ( wp ) function dnrm2 ( n , x , incx ) integer ( ip ), intent ( in ) :: n real ( wp ), intent ( in ) :: x ( * ) integer ( ip ), intent ( in ) :: incx integer ( ip ) :: ix real ( wp ) :: ssq , absxi , norm , scale if ( n < 1 . or . incx < 1 ) then norm = 0.0_wp else if ( n == 1 ) then norm = abs ( x ( 1 ) ) else scale = 0.0_wp ssq = 1.0_wp do ix = 1 , 1 + ( n - 1 ) * incx , incx if ( x ( ix ) /= 0.0_wp ) then absxi = abs ( x ( ix ) ) if ( scale < absxi ) then ssq = 1.0_wp + ssq * ( scale / absxi ) ** 2 scale = absxi else ssq = ssq + ( absxi / scale ) ** 2 end if end if end do norm = scale * sqrt ( ssq ) end if dnrm2 = norm end function dnrm2 !***************************************************************************** !> ! Scales a vector by a constant. subroutine dscal ( n , sa , x , incx ) integer ( ip ), intent ( in ) :: n !! the number of entries in the vector. real ( wp ) , intent ( in ) :: sa !! the multiplier. real ( wp ), intent ( inout ) :: x ( * ) !! the vector to be scaled. integer ( ip ), intent ( in ) :: incx !! the increment between successive entries of X. integer ( ip ) :: i , ix integer ( ip ) :: m if ( n <= 0 ) then return else if ( incx == 1 ) then m = mod ( n , 5 ) x ( 1 : m ) = sa * x ( 1 : m ) do i = m + 1 , n , 5 x ( i ) = sa * x ( i ) x ( i + 1 ) = sa * x ( i + 1 ) x ( i + 2 ) = sa * x ( i + 2 ) x ( i + 3 ) = sa * x ( i + 3 ) x ( i + 4 ) = sa * x ( i + 4 ) end do else if ( 0 <= incx ) then ix = 1 else ix = ( - n + 1 ) * incx + 1 end if do i = 1 , n x ( ix ) = sa * x ( ix ) ix = ix + incx end do end if end subroutine dscal #endif end module lsmrblas","tags":"","loc":"sourcefile/lsmrblas.f90.html"},{"title":"lsmrDataModule.F90 – LSMR","text":"Files dependent on this one sourcefile~~lsmrdatamodule.f90~~AfferentGraph sourcefile~lsmrdatamodule.f90 lsmrDataModule.F90 sourcefile~lsmrblas.f90 lsmrblas.F90 sourcefile~lsmrblas.f90->sourcefile~lsmrdatamodule.f90 sourcefile~lsmrmodule.f90 lsmrModule.f90 sourcefile~lsmrmodule.f90->sourcefile~lsmrdatamodule.f90 Help Graph Key Nodes of different colours represent the following: Graph Key Source File Source File This Page's Entity This Page's Entity Solid arrows point from a file to a file which it depends on. A file\nis dependent upon another if the latter must be compiled before the former\ncan be. Source Code !***************************************************************************** !> ! Defines `integer(ip)` and `real(wp)` kind parameters. ! !### Notes ! * JW: 1/24/2024 : The original version of this file was rewritten. module lsmrDataModule use iso_fortran_env implicit none private #ifdef REAL32 integer , parameter , public :: lsmr_wp = real32 !! real kind used by this module [4 bytes] #elif REAL64 integer , parameter , public :: lsmr_wp = real64 !! real kind used by this module [8 bytes] #elif REAL128 integer , parameter , public :: lsmr_wp = real128 !! real kind used by this module [16 bytes] #else integer , parameter , public :: lsmr_wp = real64 !! real kind used by this module [8 bytes] #endif integer , parameter , public :: lsmr_ip = int32 !! integer kind used by this module [4 bytes] end module lsmrDataModule","tags":"","loc":"sourcefile/lsmrdatamodule.f90.html"}]} \ No newline at end of file