Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
10 changes: 10 additions & 0 deletions LICENSE
Original file line number Diff line number Diff line change
Expand Up @@ -19,3 +19,13 @@ AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.


## License for src/modspainv.f90

The code in modspainv.f90 is based on code originally written by Karin Meyer
(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated
and modified code is re-licensed under the MIT license with express permission
from Karin Meyer. We are grateful for her permission to use and adapt this
work.

6 changes: 6 additions & 0 deletions README.md
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,12 @@ The brief documentation is available in the directory *doc* (see *mainpage.md*).


## Acknowledgements

The code in modspainv.f90 is based on code originally written by Karin Meyer
(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated and
modified code is re-licensed under the MIT license with express permission from
Karin Meyer. We are grateful for her permission to use and adapt this work.

This library was inspired by several sources:


Expand Down
6 changes: 5 additions & 1 deletion src/modspainv.f90
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,11 @@
!Support of sparse inverse of SPSD matrices by implementing the S. D. Kachman modifications
!(https://www.ars.usda.gov/ARSUserFiles/80420530/MTDFREML/MTDFMan.pdf ; Chapter 6)

!Rewritten for my purposes
!The code in modspainv.f90 is based on code originally written by Karin Meyer
!(didgeridoo.une.edu.au/womwiki/doku.php?id=fortran:fortran). The incorporated and
!modified code is re-licensed under the MIT license with express permission from
!Karin Meyer. We are grateful for her permission to use and adapt this work.


module modspainv
#if (_DP==0)
Expand Down
54 changes: 27 additions & 27 deletions src/modsparse_inv_int.f90
Original file line number Diff line number Diff line change
Expand Up @@ -92,40 +92,40 @@ subroutine gsfct(neqns, xlnz, lnz, xnzsub, nzsub, diag, iflag)
allocate (link(neqns), source=0)
allocate (temp(neqns), source=0.0_wp)
iflag = 0
!c --------------------------------------------
!c compute column l(*,j) for j = 1,...., neqns.
!c --------------------------------------------
!--------------------------------------------
!compute column l(*,j) for j = 1,...., neqns.
!--------------------------------------------
do j = 1, neqns
!c -------------------------------------------
!c for each column l(*,k) that affects l(*,j).
!c -------------------------------------------
!-------------------------------------------
!for each column l(*,k) that affects l(*,j).
!-------------------------------------------
diagj = 0.0d0
newk = link(j)
k = newk
do while (k /= 0)
newk = link(k)
!c ---------------------------------------
!c outer product modification of l(*,j) by
!c l(*,k) starting at first(k) of l(*,k).
!c ---------------------------------------
!---------------------------------------
!outer product modification of l(*,j) by
!l(*,k) starting at first(k) of l(*,k).
!---------------------------------------
kfirst = first(k)
ljk = lnz(kfirst)
diagj = diagj + ljk*ljk
istrt = kfirst + 1
istop = xlnz(k + 1) - 1
if (istop >= istrt) then
!c ------------------------------------------
!c before modification, update vectors first,
!c and link for future modification steps.
!c ------------------------------------------
!------------------------------------------
!before modification, update vectors first,
!and link for future modification steps.
!------------------------------------------
first(k) = istrt
i = xnzsub(k) + (kfirst - xlnz(k)) + 1
isub = nzsub(i)
link(k) = link(isub)
link(isub) = k
!c ---------------------------------------
!c the actual mod is saved in vector temp.
!c ---------------------------------------
!---------------------------------------
!the actual mod is saved in vector temp.
!---------------------------------------
do ii = istrt, istop
isub = nzsub(i)
temp(isub) = temp(isub) + lnz(ii)*ljk
Expand All @@ -134,23 +134,23 @@ subroutine gsfct(neqns, xlnz, lnz, xnzsub, nzsub, diag, iflag)
end if
k = newk
end do
!c ----------------------------------------------
!c apply the modifications accumulated in temp to
!c column l(*,j).
!c ----------------------------------------------
!----------------------------------------------
!apply the modifications accumulated in temp to
!column l(*,j).
!----------------------------------------------
diagj = diag(j) - diagj
if (diag(j) < tol) then
!c ------------------------------------------------------
!c error - zero diagonal element
!c ------------------------------------------------------
!------------------------------------------------------
!error - zero diagonal element
!------------------------------------------------------
diagj = 0.d0
diag(j) = 0.d0
iflag = iflag + 1
else
if (diagj < tol*diag(j)) then
!c ------------------------------------------------------
!c error - zero or negative square root in factorization.
!c ------------------------------------------------------
!------------------------------------------------------
!error - zero or negative square root in factorization.
!------------------------------------------------------
diagj = 0.d0
diag(j) = 0.d0
iflag = iflag + 1
Expand Down