!jInvert.f95 
!2026.04.19.1810cdt-    Complex Matrix Overwriting Inverter

!    Author- Jeffrey M. Setterholm, Lakeville,MN 55044 USA
! IP Status- Free source code (e.g.: post copyright)
!
!  Computer- "T4"/HP-800-G4-Mini/i7-8700T/IntelUHD630/Win10Pro-22H2
!             ^name  ^Mfgr.Id    ^chipset ^graphics   ^OS
!                /AbsoftProFortran 21.0.2/OpenGL+Glut3.6
!                 ^compiler               ^Fortran graphics
!Disclaimer:
! ***************************************************************
! **********   Individual cognition is always flawed,  **********
! **********         including yours and mine.         **********
! **********                 - So: -                   **********
! ********** Use this code & program at your own risk. **********
! ***************************************************************

!Table of Contents: ...use to search...
!Module jInvert_FoR                   <- The "Frame-of-Reference"          @ 31
!End Module jInvert_FoR                                                    @ 68
!Subroutine jInvert(Nin,C,iRW,iP)     <- the overwriting complex inverter  @ 68
!Subroutine jInvert_Report(N,C,iP)    <- overwriter progress reports       @237
!Subroutine jInvert_Check(N,Cout,iP)  <- quality confirmation              @270
!                                                   @: approximate Line#'s ^^^^

!    See "jInvert-f95.pdf" for a color-coded-syntax version of this file.
!------------------------------------------------------------------------------
!---------------------------------------------------------------------------7 9

Module jInvert_FoR     !The "Frame-of-Reference" for all that follows.

!use jInvert_FoR  !Complex overwriting Inverter variables            2026.04.19
!          ,only: ValMin,iRank,CdetN,iUsed,iRu,jCu,Temp   !<-any subset allowed

!The key variables of "jInvert_Demo" which reads "Invert-Input.csv":
  character(Len=40):: DataFileIn     = "jInvert-Input.txt"
  character(Len=40):: DataFileOut    = "jInvert-Output.txt"
  integer(4)       :: iRW             !Unit# for- file Reads & Writes
  integer(4)       :: iP              !         - writing DataFileOut
  character(Len=79):: YourDataHeader &!First line of DataFileIn: descriptiion
                                     = "Reading the First line failed."
  integer( 4)      :: Nin             !The size of Cin                     @ 98
  integer(4)       :: iRankin         !Rank of the inverse                 @213
  integer(4)       :: iRankOut

  complex(16),allocatable:: C(  :,:)  !Input -> overwritten -> output      @ 99

!The key variables of "jInvert()" are:
  integer( 4) :: N        !           @ 51
  real(16)    :: ValMin   !                                                @117
  real(16)    :: SNR       = -1._16   !Signal-to-Noise Ratio               @152
  integer(4)  :: iRank    !Tracks the growing rank of the inverse          @168
  complex(16) :: CdetN    !Abs(determinant product) for each iteration     @169
  integer(4) ,allocatable::iUsed(:)  !                                     @194
                                     !Overwriting:
  integer(4) ,allocatable::iRu(:)    !Tracks used pivot- rows              @175
  integer(4) ,allocatable::jCu(:)    !                 - columns           @178
  
  complex(16),allocatable::Temp(:)   !Temporary use vector of size Nin     ----
  complex(16),allocatable::Cin(:,:)  !Used by jInvert_Check()              @124                     
  complex(16),allocatable::CTest(:,:)!Used by jInvert_Check()              @288                    

!Contains
End Module jInvert_FoR
!---------------------------------------------------------------------------7 9



Subroutine jInvert(Nin,C,iRankOut,iRW,iP)
!2026.04.19.1815cdt JMS- Adapted for complex variables
!2006.05.25.0715cdt JMS- An instantiation of the "Setterholm Matrix Inverter"
!--Globals
 use jInvert_FoR &!Complex overwriting Inverter variables            2026.04.19
                 ,only: ValMin,iRank,CdetN,iUsed,iRu,jCu,Temp,Cin,Ctest,SNR
!--End Globals
 implicit none
 !--Arguments   These are arguments, not the Frame-of-Reference (FoR) variables
 !                with the same names, which jInvert_Demo uses.
 !              In the same way, you can use the FoR variables when calling
 !                jInvert().
  integer(4) :: Nin
  complex(16):: C(Nin,Nin)
  integer(4) :: iRankOut
  integer(4) :: iRW        !A Fortram unit# to use for reading or writing
  integer(4) :: iP
    !--Internals
  integer(4) :: N
  integer(4) ::i,iu,iu2,j,ju,ju2,L,iuB,juB
  complex(16)::Cmult
  complex(16)::ValMax,ValMax2
  integer(8) ::iBitFlag,jBitFlag
  integer(4) ::iAlloc
  complex(16):: Czero,Cone
   !-----                                                            !end defs
   if(iP>5) write(iP, "(/' {jInvert(Nin,C,iRankOut,iRW,iP):'&
           &,45('-'),t40,' [allocating] ',t63,' jInvert.f95 @ 97')")
   N = Nin           
   !---
   Czero = ( 0._16, 0._16) !Complex: zero
   Cone  = (+1._16, 0._16) !Complex: real 1.0

   allocate(iUsed(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: iUsed(N) allocation error. Halt.'
            iUsed  =   0
   allocate(iRu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: iRu(N)   allocation error. Halt.'
            iRu    =   0
   allocate(jCu(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: jCu(N)   allocation error. Halt.'
            jCu    =   0
   allocate(Temp(N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: Temp(N)  allocation error. Halt.'
            Temp   = Czero
   allocate(Cin(N,N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: Cin(N,N) allocation error. Halt.'
            Cin    = Czero
   allocate(CTest(N,N),stat=iAlloc)
   if(iAlloc.ne.0) stop 'jInverter: CTest(N,N) allocation error. Halt.'
            CTest  = Czero
   !---
            CdetN  = Czero
   !---
   Cin = C  !Saves C as Cin for use in analyzing Linear Dependence.
            !C will be overwritten with its inverse when C is full-rank.
            !"jInvert_Check()" aids interpreting less-than-full-rank results.

   !This is a binary file of N and Cin. Windows 10 security blocks ".bin" files
    if(iP>5) write(iP,"('  Writing binary file jCin.bim'      ,t75,' @129')")
    if(iP>5) write(iP,"('  Fortran arrays store in column-major order.',/)")
   open( unit=iRW, file="jCin.bim"   , action='write', form = 'unformatted' &
                 , access='sequential', status='replace')
     write(iRW) n
     write(iRW) C
   close(iRW)

   iUsed=0; iRu=0; jCu=0; ValMin=0._16; Temp=Czero; iRank=0
   do i=1,N    ;iRu(i)=i ;jCu(i)=i                ;end do!i
   if(iP.gt.5) call jInvert_Report(N,C,iP)                    !Initial printout

!  Progressive inversion - from the abs(largest) pivot down to the noise,
!  ... a loop:
   do L=1,N ;iU=0     ;jU=0     ;ValMax =Czero  
     do i=1,N   ;if(iRu(i).lt.0)                                          cycle
       do j=1,N ;if(jCu(j).lt.0)                                          cycle
         if(abs(ValMin).gt.abs(C(i,j)))                                   cycle
         if(abs(ValMax).gt.abs(C(i,j)))                                   cycle
                ValMax=C(i,j) ;iu=i ;ju=j
       enddo!j
     enddo!i
     if(L==1) then
       ValMin=abs(ValMax)*1.e-24_16 !"noise floor" ...not cast in concrete.
                                  ! Works here because Fortran's complex(16)
                                  ! components have 31-decimal-place precision.
                                  !Can be overwritten by SNR...
       if(SNR>0._16) then
         ValMin=abs(ValMax)/SNR                                     !2026.02.24
         if(iP>5) &
         write(iP,"(34x,'User SNR =',e15.6,' Signal-to-Noise Ratio')") SNR
       endif!(SNR>0._16)
       if(iP>5) &
         write(iP,"(34x,' ValMin  =',e15.6,' scaled by pivot:'/)") ValMin
       CdetN=ValMax
     endif!(L==1)  

     if(iU==0)                                                             Exit
     if(abs(ValMax).lt.ValMin)                                             Exit
     iRank=iRank+1
     if(L>1) CdetN=CdetN*ValMax
     ValMax=C(iu,ju) ;iu2=jCu(ju) ;ju2=iRu(iu)
     if(iP.gt.5) write(iP, &
    "('    iu =',i3,'  ju = ',i3,16x,'Pivot =',sp,' (',f13.6,',',f13.6,'j)')")&
               iu      ,  ju,                   ValMax
     Temp(1:N)=C(1:N,ju)                            ;i=iRu(iu)
               C(1:N,ju)=C(1:N,ju2)                 ;  iRu(iu)=iRu(iu2)
                         C(1:N,ju2)=Temp(1:N)       ;          iRu(iu2)=-abs(i)
     Temp(1:N)=C(iu,1:N)                            ;j=jCu(ju)
               C(iu,1:N)=C(iu2,1:N)                 ;  jCu(ju)=jCu(ju2)
                         C(iu2,1:N)=Temp(1:N)/ValMax ;         jCu(ju2)=-abs(j)
     do i=1,N ;if(i.eq.iu2) cycle ;Cmult=C(i,ju2)
       do j=1,N ;C(i,j)=C(i,j)-C(iu2,j)*Cmult ;enddo!j
       C(i,ju2)=-Cmult/ValMax
     enddo; C(iu2,ju2)=1._16/ValMax
     if(iP.gt.5) call jInvert_Report(N,C,iP)                     !Loop printout
   enddo !Look for the next largest remaining pivot.
   
  !Done trying to invert C. Report how it turned out:
   if(iP>5) write(iP, &
     "(13x,'sign untracked:',9x,'CdetN =',sp,' (',f13.6,',',f13.6,'j)')") CdetN
                                      iBitFlag=0; jBitFlag=0
   do i=1,N
      if((iRu(i).lt.0).and.(i.lt.62)) iBitFlag=iBitFlag+2**(i-1)
      if((jCu(i).lt.0).and.(i.lt.62)) jBitFlag=jBitFlag+2**(i-1)
      if( iRu(i).lt.0)iUsed(i)=1
!     Zero linearly-dependent rows and columns, if any:
      if(iRu(i).gt.0) C(i,1:N)=Czero
      if(jCu(i).gt.0) C(1:N,i)=Czero
   enddo!i

   if((iP >5).and.(iRank.lt.N)) then
     write(iP," (/'  ^:Linearly Dependent. Rows and Columns with&
               & positive indices (above) failed to invert,')")
     write(iP," (' The corresponding rows & columns are zeroed.',\)")
     write(iP," (' Bitflags: 1-used, 0-discarded [N->1]')")        
     write(iP,"('[',b63.6,'] input' )") iBitFlag !For first 61
     write(iP,"('[',b63.6,'] output')") jBitFlag !For first 61
     call jInvert_Report(N,C,iP)  !Linearly-dependence printout
   endif!iP>5 & iu==0

  !if(iRank<Nin) call jInvert_Check(N,C,iP)
                 call jInvert_Check(N,C,iP)
   
    iRankOut = iRank
   !This is a binary file of N iRank, and Cout.
    if(iP>5) write(iP,"(/'  Writing binary file jCinv.bim '     ,t75,' @214')")
    if(iP>5) write(iP,"( '  Fortran arrays store in column-major order.')")
    open( unit=iRW, file="jCinv.bim"   , action='write', form = 'unformatted' &
                  , access='sequential', status='replace')
     write(iRW) n
     write(iRW) iRank
     write(iRW) C
   close(iRW)

  !Housekeeping:
   deallocate(iUsed)
   deallocate(iRu)
   deallocate(jCu)
   deallocate(Temp)
   deallocate(Cin)
   deallocate(CTest)
    if(iP>5) write(iP,&
    "('  jInvert} ',60('-'),t40,' [deallocated] '        ,t70,' Exit @230',/)")
                                                                         return
End Subroutine jInvert
!---------------------------------------------------------------------------7-9

Subroutine jInvert_Report(N,C,iP)
!2026.01.31.0745cst JMS- For complex matrices
!2018.09.04.1340cdt JMS- Prints Invert's progress.
!                          Expects,uses, & returns real(16)'s (~QUAD precision).
 !---                                                                No globals
!--Globals
 use jInvert_FoR &!Complex overwriting Inverter variables            2026.04.19
                 ,only: ValMin,iRank,CdetN,iRu,jCu
!--End Globals
 implicit none
 !--Arguments
  integer(4) :: N
  complex(16):: C(N,N)
  real(16)   :: Noise
  integer(4) :: iP
 !--Internals
  integer(4)::i,j
  !--EndDefs----------------------------------------------
   Noise = ValMin                                            ;if(iP < 6) return
   write(iP,"('@ Rank =',i3,21x,'abs(CdetN) =',e23.10/22x,' ',\)") &
               iRank,abs(CdetN)
   do j=1,N;   write(iP,"(sp,i6,24x,' ',\)") jCu(j)
   enddo;      write(iP,"('')")
   do i=1,N;   write(iP,"(i,' ' ,\)") iRu(i)
     do j=1,N; write(iP,"(sp,' (',f13.6,',',f13.6,'j)',\)") C(i,j)
     enddo;    write(iP,"('')")
   enddo!i
             write(iP,"(13(' '),\)")  
   do i=1,N; write(iP,"(31('-'),\)")
   enddo;    write(iP,*);                                                return
End Subroutine jInvert_Report
!---------------------------------------------------------------------------7-9

Subroutine jInvert_Check(N,Cout,iP)
!2026.04.19.1605cdt JMS- Prints C*Cin and Cin*C
!--Globals
 use jInvert_FoR &!Complex overwriting Inverter variables            2026.04.19
                 ,only: Cin,CTest,iRank
!--End Globals
 implicit none
 !--Arguments
  integer(4) :: N
  complex(16):: Cout(N,N)
  integer(4) :: iP
 !--Internals
  complex(16):: Cz = ( 0., 0.)
  integer(4)::i,j,k
  !--EndDefs----------------------------------------------
   CTest = Cz
   do i = 1,N
     do j = 1,N
       do k = 1,N
         CTest(i,j) = CTest(i,j) + Cout(i,k)*Cin(k,j)
       enddo!k
     enddo!j
   enddo!i
   if(iP>5) write(iP,"(/'Cinv*Cin: maps inputs back into inputs.'&
                       &,' An Identity matrix if full-rank.')")
   do i=1,N
     do j=1,N; write(iP,"(sp,' (',e13.4,',',e13.4,'j),',\)") CTest(i,j)
     enddo;    write(iP,"('')")
   enddo!i

   CTest = Cz
   do i = 1,N
     do j = 1,N
       do k = 1,N
         CTest(i,j) = CTest(i,j) + Cin(i,k)*Cout(k,j)
       enddo!k
     enddo!j
   enddo!i
   if(iP>5) write(iP,"(/'Cin*Cinv: maps outputs back into outputs.'&
                       &,' An Identity matrix if full-rank.')")
                       
                       
   do i=1,N
     do j=1,N; write(iP,"(sp,' (',e13.4,',',e13.4,'j),',\)") CTest(i,j)
     enddo;    write(iP,"('')")
   enddo!i
   if((iP>5).and.(iRank<N)) &
   write(iP,"( /'Consider these results tentative for linearly dependent&
              & complex matrices.',/,'  I have not applied them yet.'&
              &,/'For real matrices the results have confirmed utility.')")
                                                                         return
End Subroutine jInvert_Check
!---------------------------------------------------------------------------7-9
