diff --git a/CMakeLists.txt b/CMakeLists.txt
index fe41418..7d72a97 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -577,8 +577,14 @@ if(EXT_MODTEST)
   add_target_exe_serial_wrapper(test_zgetrf test_tools)
   add_target_exe_serial_wrapper(test_dgetrfi test_tools)
   add_target_exe_serial_wrapper(test_zgetrfi test_tools)
+  add_target_exe_serial_wrapper(test_zsktri test_tools)
+  add_target_exe_serial_wrapper(test_zsktrs test_tools)
   add_target_exe_serial_wrapper(test_openmp_reduction_real test_tools)
   add_target_exe_serial_wrapper(test_openmp_reduction_complex test_tools)
+  add_target_exe_serial_wrapper(test_upwinvp test_tools)
+  add_target_exe_serial_wrapper(test_upwinvp_complex test_tools)
+  #add_target_exe_serial_wrapper(test_upwinvp_pfaff test_tools)
+  #add_target_exe_serial_wrapper(test_upwinvp_pfaff_complex test_tools)
 
   if(EXT_GPU)
 
diff --git a/src/a_module_tests/CMakeLists.txt b/src/a_module_tests/CMakeLists.txt
index 79ce57f..8c3294d 100644
--- a/src/a_module_tests/CMakeLists.txt
+++ b/src/a_module_tests/CMakeLists.txt
@@ -162,6 +162,44 @@ foreach( EXECUTABLE IN LISTS EXECUTABLES_S_L
       )
   endif()
 
+  ######################################################################
+  #
+  # Add sources and link libraries for zsktr[ri] wrapper tests
+  #
+
+  if(${EXECUTABLE} MATCHES test_zsktr[si])
+      set_target_properties( ${EXECUTABLE} PROPERTIES SUFFIX ".x")
+    
+    target_sources( ${EXECUTABLE}
+      PRIVATE
+      ${EXECUTABLE}.f90
+      )
+    
+    target_link_libraries( ${EXECUTABLE}
+      PRIVATE
+        common-serial
+      )
+  endif()
+
+  ######################################################################
+  #
+  # Add sources and link libraries for upwinvp wrapper tests
+  #
+
+  if(${EXECUTABLE} MATCHES test_upwinvp.*)
+      set_target_properties( ${EXECUTABLE} PROPERTIES SUFFIX ".x")
+    
+    target_sources( ${EXECUTABLE}
+      PRIVATE
+      ${EXECUTABLE}.f90
+      )
+    
+    target_link_libraries( ${EXECUTABLE}
+      PRIVATE
+        common-serial
+      )
+  endif()
+
   ######################################################################
   #
   #
diff --git a/src/a_module_tests/test_upwinvp.f90 b/src/a_module_tests/test_upwinvp.f90
new file mode 100644
index 0000000..d9de578
--- /dev/null
+++ b/src/a_module_tests/test_upwinvp.f90
@@ -0,0 +1,53 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+program test_upwinvp
+    use constants, only: yes_ontarget
+    implicit none
+    integer, parameter :: nel = 3, indt = 4
+    real*8 :: psi(indt, nel, 2), ainv(nel), ainvn(nel), winv(nel, indt)
+    integer :: i, j, k
+    real*8, parameter :: eps = 1.0e-10 ! A small value for numerical comparison
+
+    yes_ontarget = .true.
+
+    ! Initialize test data
+    winv = 0.0
+    do i = 1, indt
+        do j = 1, nel
+            do k = 1, 2
+                psi(i, j, k) = i*j*k
+            end do
+        end do
+    end do
+    ainv = (/1.0, 2.0, 3.0/)
+    ainvn = (/4.0, 5.0, 6.0/)
+
+    call upwinvp(nel, indt, winv, ainv, ainvn, psi)
+
+    ! Check the result
+    do i = 1, indt
+        do j = 1, nel
+            if (abs(winv(j, i) - (ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2))) > eps) then
+                print *, "Test failed at i=", i, ", j=", j
+                print *, "ERROR"
+                stop 1
+            end if
+        end do
+    end do
+
+    print *, "OK"
+    print *, "All tests passed."
+end program test_upwinvp
diff --git a/src/a_module_tests/test_upwinvp_complex.f90 b/src/a_module_tests/test_upwinvp_complex.f90
new file mode 100644
index 0000000..766ac02
--- /dev/null
+++ b/src/a_module_tests/test_upwinvp_complex.f90
@@ -0,0 +1,60 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+program test_upwinvp_complex
+    use constants, only: yes_ontarget
+    implicit none
+    integer, parameter :: nel = 3, indt = 4
+    complex*16 :: psi(indt, nel, 2), ainv(nel), ainvn(nel), winv(nel, indt), expected(nel, indt)
+    integer :: i, j, k
+    real*8, parameter :: eps = 1.0e-10 ! A small value for numerical comparison
+
+    yes_ontarget = .true.
+
+    ! Initialize test data
+    winv = cmplx(0.0, 0.0)
+    do i = 1, indt
+        do j = 1, nel
+            do k = 1, 2
+                psi(i, j, k) = cmplx(i*j*k, i*j*k)
+            end do
+        end do
+    end do
+    ainv = (/cmplx(1.0, 1.0), cmplx(2.0, 2.0), cmplx(3.0, 3.0)/)
+    ainvn = (/cmplx(4.0, 4.0), cmplx(5.0, 5.0), cmplx(6.0, 6.0)/)
+
+    call upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
+
+    ! Calculate expected result
+    do i = 1, indt
+        do j = 1, nel
+            expected(j, i) = ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
+        end do
+    end do
+
+    ! Check the result
+    do i = 1, indt
+        do j = 1, nel
+            if (abs(winv(j, i) - expected(j, i)) > eps) then
+                print *, "Test failed at i=", i, ", j=", j
+                print *, "OK"
+                stop 1
+            end if
+        end do
+    end do
+
+    print *, "OK"
+    print *, "All tests passed."
+end program test_upwinvp_complex
diff --git a/src/a_module_tests/test_zsktri.f90 b/src/a_module_tests/test_zsktri.f90
new file mode 100644
index 0000000..5df78b9
--- /dev/null
+++ b/src/a_module_tests/test_zsktri.f90
@@ -0,0 +1,96 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+program test_zsktri
+    implicit none
+
+    complex*16, allocatable, dimension(:, :) :: A, A_inv, A_inv_orig
+    complex*16, allocatable, dimension(:) :: W
+    integer, allocatable, dimension(:) :: ipiv
+    real*8, allocatable, dimension(:, :) :: helper_r, helper_c
+    complex*16 :: one = 1.d0, zero = 0.d0
+    integer :: s, gen, ii, jj, info
+    character(len=1) :: uplo
+
+    ! s dimension of the test matrix, s x s.
+    ! uplo: U->upper triangular, L->lower triangular
+    ! gen = 0 : Compare matrices, gen = 1 : Generate matrices
+    write (*, *) 's, uplo, gen'
+    read (*, *) s, uplo, gen
+
+    allocate (ipiv(s))
+    allocate (A(s, s))
+    allocate (A_inv(s, s))
+    allocate (A_inv_orig(s, s))
+    allocate (W(s**2 + 12*s - 2))
+    allocate (helper_r(s, s))
+    allocate (helper_c(s, s))
+
+    if (gen .eq. 1) then
+
+        ! generate a skew_symmetric matrix A
+        call random_number(helper_r)
+        call random_number(helper_c)
+
+        do ii = 1, s
+            do jj = 1, ii - 1
+                helper_c(jj, ii) = -helper_c(ii, jj)
+                helper_r(jj, ii) = -helper_r(ii, jj)
+            end do
+        end do
+
+        do ii = 1, s
+            helper_r(ii, ii) = 0
+            helper_c(ii, ii) = 0
+        end do
+
+        A = cmplx(helper_r, helper_c)
+
+    else
+        open (unit=10, form="unformatted", file="A", action="read")
+        read (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="A_inv_orig", action="read")
+        read (10) A_inv_orig
+        close (10)
+    end if
+
+    W = 0
+
+    do ii = 1, s
+        ipiv(ii) = ii
+    end do
+
+    call zsktri(uplo, s, A, s, A_inv, s, ipiv, W, info)
+
+    if (gen .eq. 1) then
+        open (unit=10, form="unformatted", file="A", action="write")
+        write (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="A_inv_orig", action="write")
+        write (10) A_inv_orig
+        close (10)
+    else
+        A_inv = A_inv - A_inv_orig
+        if (maxval(abs(A_inv)) > 1.0d-10) then
+            print *, "ERROR"
+        else
+            print *, "OK"
+        end if
+    end if
+
+    deallocate (A, A_inv, A_inv_orig, W, helper_r, helper_c)
+
+end program test_zsktri
diff --git a/src/a_module_tests/test_zsktrs.f90 b/src/a_module_tests/test_zsktrs.f90
new file mode 100644
index 0000000..7c058c9
--- /dev/null
+++ b/src/a_module_tests/test_zsktrs.f90
@@ -0,0 +1,77 @@
+! Copyright (C) 2022 TurboRVB group
+!
+! This program is free software: you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation, either version 3 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License
+! along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+program test_zsktrs
+    implicit none
+
+    complex*16, allocatable, dimension(:, :) :: B, B_orig, X
+    complex*16, allocatable, dimension(:) :: A
+    real*8, allocatable, dimension(:) :: helper_r, helper_c
+    complex*16 :: one = 1.d0, zero = 0.d0
+    integer :: s, gen, ii, jj, info
+    character(len=1) :: uplo
+
+    ! s dimension of the test matrix, s x s.
+    ! uplo: U->upper triangular, L->lower triangular
+    ! gen = 0 : Compare matrices, gen = 1 : Generate matrices
+    write (*, *) 's, uplo, gen'
+    read (*, *) s, uplo, gen
+
+    allocate (A(s - 1))
+    allocate (B(s, s))
+    allocate (B_orig(s, s))
+    allocate (X(s, s))
+    allocate (helper_r(s - 1))
+    allocate (helper_c(s - 1))
+
+    if (gen .eq. 1) then
+
+        ! generate tridiagonal skewsymmetric matrix A
+        call random_number(helper_r)
+        call random_number(helper_c)
+        A = cmplx(helper_r, helper_c)
+
+    else
+        open (unit=10, form="unformatted", file="A", action="read")
+        read (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="B_orig", action="read")
+        read (10) B_orig
+        close (10)
+    end if
+
+    X = 0
+
+    call zsktrs(uplo, s, s, A, B, s, X, info)
+
+    if (gen .eq. 1) then
+        open (unit=10, form="unformatted", file="A", action="write")
+        write (10) A
+        close (10)
+        open (unit=10, form="unformatted", file="B_orig", action="write")
+        write (10) B_orig
+        close (10)
+    else
+        B = B - B_orig
+        if (maxval(abs(B)) > 1.0d-10) then
+            print *, "ERROR"
+        else
+            print *, "OK"
+        end if
+    end if
+
+    deallocate (A, B, B_orig, X, helper_r, helper_c)
+
+end program test_zsktrs
diff --git a/src/m_common/upwinv.f90 b/src/m_common/upwinv.f90
index b3abd1f..457529f 100644
--- a/src/m_common/upwinv.f90
+++ b/src/m_common/upwinv.f90
@@ -16,8 +16,15 @@
 subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    real*8 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    real*8, intent(in) :: psi(indt, nel), v(nel)
+    real*8, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
+
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif
@@ -43,11 +50,18 @@ subroutine upwinv(nel, jel, indt, nelorb, winv, v, psi)
     end do
     return
 end
+
 subroutine upwinv_complex(nel, jel, indt, nelorb, winv, v, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j, jel, nelorb
-    complex*16 winv(nel, indt), psi(indt, nel), v(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt, jel, nelorb
+    complex*16, intent(in) :: psi(indt, nel), v(nel)
+    complex*16, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
 #ifdef _OFFLOAD
 !$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
 #endif
diff --git a/src/m_common/upwinvp.f90 b/src/m_common/upwinvp.f90
index 3062009..fe835f2 100644
--- a/src/m_common/upwinvp.f90
+++ b/src/m_common/upwinvp.f90
@@ -13,11 +13,19 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> subroutine to update the matrix winv, for real(8)
 subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
     use constants, only: yes_ontarget
     implicit none
-    integer nel, indt, i, j
-    real*8 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt
+    real*8, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
+    real*8, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer :: i, j
+
 #ifdef _OFFLOAD
     if (yes_ontarget) then
 !$omp target teams distribute parallel do collapse(2)
@@ -42,12 +50,48 @@ subroutine upwinvp(nel, indt, winv, ainv, ainvn, psi)
 #endif
     return
 end
+
+!> subroutine to update the matrix winv, for complex(16)
+subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
+    use constants, only: yes_ontarget
+    implicit none
+
+    ! argument parameters
+    integer, intent(in) :: nel, indt
+    complex*16, intent(in) :: psi(indt, nel, 2), ainv(nel), ainvn(nel)
+    complex*16, intent(inout) :: winv(nel, indt)
+
+    ! local variables
+    integer i, j
+
+#ifdef _OFFLOAD
+!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
+#endif
+    do i = 1, indt
+        do j = 1, nel
+            winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
+        end do
+    end do
+#ifdef _OFFLOAD
+!$omp end target teams distribute parallel do
+#endif
+    return
+end
+
+!> subroutine to update the matrix winv of Pfaffian, for real(8)
 subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
     use constants, only: yes_ontarget
     implicit none
-    integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
-    real*8 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)
+
+    ! argument parameters
+    integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
+    real*8, intent(in) :: psi(nmolipf, indt), ainv(nmol)
+    real*8, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)
+
+    ! local variables
+    integer i, j, nelcdo
     real*8 csum
+
     if (yes_ontarget) then
     if (nelc .le. nelup) then
 #ifdef _OFFLOAD
@@ -116,11 +160,19 @@ subroutine upwinvp_pfaff(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, win
     end if
     return
 end
+
+!> subroutine to update the matrix winv of Pfaffian, for complex(16)
 subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt, winvup, winvdo, psi, ainv)
     use constants, only: yes_ontarget
     implicit none
-    integer nelc, nelcdo, nelup, neldo, nmol, nmolipf, nmolshift, indt, i, j
-    complex*16 winvup(nelup, indt), winvdo(neldo, indt), psi(nmolipf, indt), ainv(nmol)
+
+    ! argument parameters
+    integer, intent(in) :: nelc, nelup, neldo, nmol, nmolipf, nmolshift, indt
+    complex*16, intent(in) :: psi(nmolipf, indt), ainv(nmol)
+    complex*16, intent(inout) :: winvup(nelup, indt), winvdo(neldo, indt)
+
+    ! local variables
+    integer i, j, nelcdo
     complex*16 csum
     if (yes_ontarget) then
         if (nelc .le. nelup) then
@@ -194,21 +246,3 @@ subroutine upwinvp_pfaff_complex(nelc, nelup, neldo, nmol, nmolipf, nmolshift, i
     end if
     return
 end
-subroutine upwinvp_complex(nel, indt, winv, ainv, ainvn, psi)
-    use constants, only: yes_ontarget
-    implicit none
-    integer nel, indt, i, j
-    complex*16 winv(nel, indt), psi(indt, nel, 2), ainv(nel), ainvn(nel)
-#ifdef _OFFLOAD
-!$omp target teams distribute parallel do collapse(2) if(yes_ontarget)
-#endif
-    do i = 1, indt
-        do j = 1, nel
-            winv(j, i) = winv(j, i) + ainv(j)*psi(i, j, 1) + ainvn(j)*psi(i, j, 2)
-        end do
-    end do
-#ifdef _OFFLOAD
-!$omp end target teams distribute parallel do
-#endif
-    return
-end
diff --git a/src/m_common/write_type_orb.f90 b/src/m_common/write_type_orb.f90
index 1702557..0fbd79f 100644
--- a/src/m_common/write_type_orb.f90
+++ b/src/m_common/write_type_orb.f90
@@ -13,9 +13,15 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine returns typeorb for jastrow orbitals
 subroutine write_type_orb(nshellj, multij, ioccj, typeorb)
     implicit none
-    integer nshellj, multij(*), ioccj(*), typeorb(*)
+
+    ! argument parameters
+    integer, intent(in) :: nshellj, multij(*), ioccj(*)
+    integer, intent(out) :: typeorb(*)
+
+    ! local variables
     integer i, ii, j, ind_type
 
     ii = 0
diff --git a/src/m_common/zgemm_my.f90 b/src/m_common/zgemm_my.f90
index ba3baa0..571c505 100644
--- a/src/m_common/zgemm_my.f90
+++ b/src/m_common/zgemm_my.f90
@@ -13,6 +13,7 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine is a wrapper for the LAPACK zgemm routine
 subroutine ZGEMM_MY(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA&
         &, B, LDB, BETA, C, LDC, nproc, rank, comm_mpi)
     implicit none
diff --git a/src/m_common/zsktri.f90 b/src/m_common/zsktri.f90
index 4281420..855b12d 100644
--- a/src/m_common/zsktri.f90
+++ b/src/m_common/zsktri.f90
@@ -13,19 +13,29 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine manipulates a skewsymmetric matrix A.
 subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
     implicit none
-    character uplo
-    integer n, i, j, lda, ldinv, info, lwork, ipiv(*)
-    complex*16 a(lda, *), ainv(ldinv, *), work(*)
+
+    ! argument parameters
+    character, intent(in) :: uplo
+    integer, intent(in) :: n, lda, ldinv, ipiv(*)
+    complex*16, intent(in out) :: a(lda, *), work(*)
+    complex*16, intent(out) :: ainv(ldinv, *)
+    integer, intent(out) :: info
+
+    ! local variables
+    integer i, j
     real*8 rcond
     logical yeslap
+
     !     First factorization assumed
-    !      lwork=3*n
+    !     lwork=3*n
     !     CALL ZSKTRF( UPLO, 'N', N, A, LDA, IPIV, WORK, LWORK, INFO)
     !     Identity matrix
     !     ipiv dimension required 3n
     !     work dimension required n^2+12*n-2
+    info = 0
     yeslap = .false. ! if false the homemade algorithm is done.
     do i = 1, n
         do j = 1, i - 1
@@ -36,13 +46,13 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
             ainv(j, i) = dcmplx(0.d0, 0.d0)
         end do
     end do
-    !     fisrt permutation of ainv
+    ! fisrt permutation of ainv
     do i = n, 1, -1
         work(1:n) = ainv(i, 1:n)
         ainv(i, 1:n) = ainv(ipiv(i), 1:n)
         ainv(ipiv(i), 1:n) = work(1:n)
     end do
-    !      SECOND
+    ! second
     if (UPLO .eq. 'u' .or. UPLO .eq. 'U') then
         do i = 1, N - 1
             work(i) = a(i, i + 1)
@@ -63,7 +73,7 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
         a(2:n, 1) = dcmplx(0.d0, 0.d0)
     end if
     work(2*N - 1:3*N - 2) = dcmplx(0.d0, 0.d0) ! diagonal elements of skew matrix , obviously set to zero.
-    call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO)
+    call ZTRTRS(UPLO, 'N', 'U', N, N, A, LDA, ainv, LDINV, INFO) ! ainv = A^-1
     ! USE STANDARD LAPACK DGTSVX or the simples DGTSV?
     if (yeslap) then
         call ZGTSVX('N', 'N', N, N, work(N), work(2*N - 1), work, work(3*N)&
@@ -83,7 +93,6 @@ subroutine zsktri(uplo, n, a, lda, ainv, ldinv, ipiv, work, info)
         ainv(i, 1:n) = ainv(ipiv(i), 1:n)
         ainv(ipiv(i), 1:n) = work(1:n)
     end do
-
     return
 end
 
diff --git a/src/m_common/zsktrs.f90 b/src/m_common/zsktrs.f90
index 2effe03..638a85e 100644
--- a/src/m_common/zsktrs.f90
+++ b/src/m_common/zsktrs.f90
@@ -1,4 +1,4 @@
-! Copyright (C) 2022 TurboRVB group
+! Copyright (C) 2023- TurboRVB group
 !
 ! This program is free software: you can redistribute it and/or modify
 ! it under the terms of the GNU General Public License as published by
@@ -13,11 +13,22 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine solves A X = B
+!> To be deleted!
 subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
     implicit none
-    integer n, nhrs, ldb, info, i, j
-    complex*16 a(*), b(ldb, nhrs), x(nhrs, *)
-    character uplo
+
+    ! argument parameters
+    integer, intent(in) :: n, nhrs, ldb
+    integer, intent(out) :: info
+    complex*16, intent(in out) :: a(*), b(ldb, nhrs)
+    complex*16, intent(out) :: x(nhrs, *)
+    character, intent(in) :: uplo
+
+    ! local variables
+    integer i, j
+
+    info = 0
     if (uplo .eq. 'l' .or. uplo .eq. 'L') then
         a(1:n - 1) = -a(1:n - 1)
     end if
diff --git a/src/m_common/zsktrs_qp.f90 b/src/m_common/zsktrs_qp.f90
index aa0c9f9..272d47a 100644
--- a/src/m_common/zsktrs_qp.f90
+++ b/src/m_common/zsktrs_qp.f90
@@ -13,16 +13,31 @@
 ! You should have received a copy of the GNU General Public License
 ! along with this program. If not, see <http://www.gnu.org/licenses/>.
 
+!> This subroutine solves A * X = B.
+!> Our homemade algorithm to solve A * X = B (No transpose)
+!> where A is a skewsymmetric matrix, X and B are rectanglar matrices.
+!> The input vector a stores the elements of A such that
+!> A(i+1,i)=a(i) and A(i,i+1)=-a(i) for i=1,..,n-1
+!> a(n), b(ldb, n), x(n)
 subroutine zsktrs(uplo, n, nhrs, a, b, ldb, x, info)
     implicit none
-    integer n, nhrs, ldb, info, i, j
-    complex*16 a(*), b(ldb, *), x(*)
+
+    ! argument variables
+    integer, intent(in) :: n, nhrs, ldb
+    integer, intent(out) :: info
+    complex*16, intent(inout) :: a(*), b(ldb, *)
+    complex*16, intent(out) :: x(*)
+    character, intent(in) :: uplo
+
+    ! local variables
+    integer i, j
 #ifdef __PORT
     complex*16 aq, a2q, xo
 #else
     complex*32 aq, a2q, xo
 #endif
-    character uplo
+
+    info = 0
     if (uplo .eq. 'l' .or. uplo .eq. 'L') then
         a(1:n - 1) = -a(1:n - 1)
     end if
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index eee8aa6..610300b 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -129,7 +129,18 @@ if(EXT_MODTEST)
   add_subdirectory(test_zgetrf_128)
   add_subdirectory(test_dgetrfi_128)
   add_subdirectory(test_zgetrfi_128)
+  
+  add_subdirectory(test_zsktri_128_L)
+  add_subdirectory(test_zsktri_128_U)
+
+  add_subdirectory(test_zsktrs_128_L)
+  add_subdirectory(test_zsktrs_128_U)
 
+  add_subdirectory(test_upwinvp)
+  add_subdirectory(test_upwinvp_complex)
+  #add_subdirectory(test_upwinvp_pfaff)
+  #add_subdirectory(test_upwinvp_pfaff_complex)
+  
   #
   ######################################################################
 
diff --git a/test/test_upwinvp/CMakeLists.txt b/test/test_upwinvp/CMakeLists.txt
new file mode 100644
index 0000000..7d6f966
--- /dev/null
+++ b/test/test_upwinvp/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test upwinvp wrapper"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_upwinvp> test_upwinvp.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_upwinvp/cm.test.sh b/test/test_upwinvp/cm.test.sh
new file mode 100755
index 0000000..7676287
--- /dev/null
+++ b/test/test_upwinvp/cm.test.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	UPWINVP=$1
+	OUT=$2
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$UPWINVP" ]; then
+    echo "Executable $UPWINVP does not exists"
+    exit 1
+fi
+
+echo " Test upwinvp wrapper " 
+echo 1 | $UPWINVP > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_upwinvp_complex/CMakeLists.txt b/test/test_upwinvp_complex/CMakeLists.txt
new file mode 100644
index 0000000..a1cf5a8
--- /dev/null
+++ b/test/test_upwinvp_complex/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test upwinvp_complex wrapper"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_upwinvp_complex> test_upwinvp_complex.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_upwinvp_complex/cm.test.sh b/test/test_upwinvp_complex/cm.test.sh
new file mode 100755
index 0000000..a4620e1
--- /dev/null
+++ b/test/test_upwinvp_complex/cm.test.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	UPWINVP_COMPLEX=$1
+	OUT=$2
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$UPWINVP_COMPLEX" ]; then
+    echo "Executable $UPWINVP_COMPLEX does not exists"
+    exit 1
+fi
+
+echo " Test upwinvp_complex wrapper " 
+echo 1 | $UPWINVP_COMPLEX > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_zsktri_128_L/A b/test/test_zsktri_128_L/A
new file mode 100644
index 0000000..a2e66ce
Binary files /dev/null and b/test/test_zsktri_128_L/A differ
diff --git a/test/test_zsktri_128_L/A_inv_orig b/test/test_zsktri_128_L/A_inv_orig
new file mode 100644
index 0000000..eb44548
Binary files /dev/null and b/test/test_zsktri_128_L/A_inv_orig differ
diff --git a/test/test_zsktri_128_L/CMakeLists.txt b/test/test_zsktri_128_L/CMakeLists.txt
new file mode 100644
index 0000000..82f4c19
--- /dev/null
+++ b/test/test_zsktri_128_L/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test zsktri wrapper L (128)"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_zsktri> L test_zsktri.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_zsktri_128_L/cm.test.sh b/test/test_zsktri_128_L/cm.test.sh
new file mode 100755
index 0000000..9b3f47b
--- /dev/null
+++ b/test/test_zsktri_128_L/cm.test.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	ZSKTRI=$1
+	UPLO=$2
+	OUT=$3
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$ZSKTRI" ]; then
+    echo "Executable $ZSKTRI does not exists"
+    exit 1
+fi
+
+echo " Test zsktri wrapper " 
+cat s$UPLO | $ZSKTRI > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_zsktri_128_L/sL b/test/test_zsktri_128_L/sL
new file mode 100644
index 0000000..463f8b5
--- /dev/null
+++ b/test/test_zsktri_128_L/sL
@@ -0,0 +1 @@
+128 L 0
diff --git a/test/test_zsktri_128_U/A b/test/test_zsktri_128_U/A
new file mode 100644
index 0000000..51a66cd
Binary files /dev/null and b/test/test_zsktri_128_U/A differ
diff --git a/test/test_zsktri_128_U/A_inv_orig b/test/test_zsktri_128_U/A_inv_orig
new file mode 100644
index 0000000..eb44548
Binary files /dev/null and b/test/test_zsktri_128_U/A_inv_orig differ
diff --git a/test/test_zsktri_128_U/CMakeLists.txt b/test/test_zsktri_128_U/CMakeLists.txt
new file mode 100644
index 0000000..c143b68
--- /dev/null
+++ b/test/test_zsktri_128_U/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test zsktri wrapper U (128)"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_zsktri> U test_zsktri.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_zsktri_128_U/cm.test.sh b/test/test_zsktri_128_U/cm.test.sh
new file mode 100755
index 0000000..9b3f47b
--- /dev/null
+++ b/test/test_zsktri_128_U/cm.test.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	ZSKTRI=$1
+	UPLO=$2
+	OUT=$3
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$ZSKTRI" ]; then
+    echo "Executable $ZSKTRI does not exists"
+    exit 1
+fi
+
+echo " Test zsktri wrapper " 
+cat s$UPLO | $ZSKTRI > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_zsktri_128_U/sU b/test/test_zsktri_128_U/sU
new file mode 100644
index 0000000..95cd965
--- /dev/null
+++ b/test/test_zsktri_128_U/sU
@@ -0,0 +1 @@
+128 U 0
diff --git a/test/test_zsktrs_128_L/A b/test/test_zsktrs_128_L/A
new file mode 100644
index 0000000..077e405
Binary files /dev/null and b/test/test_zsktrs_128_L/A differ
diff --git a/test/test_zsktrs_128_L/B_orig b/test/test_zsktrs_128_L/B_orig
new file mode 100644
index 0000000..eb44548
Binary files /dev/null and b/test/test_zsktrs_128_L/B_orig differ
diff --git a/test/test_zsktrs_128_L/CMakeLists.txt b/test/test_zsktrs_128_L/CMakeLists.txt
new file mode 100644
index 0000000..3e5f50e
--- /dev/null
+++ b/test/test_zsktrs_128_L/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test zsktrs wrapper L (128)"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_zsktrs> L test_zsktrs.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_zsktrs_128_L/cm.test.sh b/test/test_zsktrs_128_L/cm.test.sh
new file mode 100755
index 0000000..f72b445
--- /dev/null
+++ b/test/test_zsktrs_128_L/cm.test.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	ZSKTRS=$1
+	UPLO=$2
+	OUT=$3
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$ZSKTRS" ]; then
+    echo "Executable $ZSKTRS does not exists"
+    exit 1
+fi
+
+echo " Test zsktrs wrapper " 
+cat s$UPLO | $ZSKTRS > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_zsktrs_128_L/sL b/test/test_zsktrs_128_L/sL
new file mode 100644
index 0000000..463f8b5
--- /dev/null
+++ b/test/test_zsktrs_128_L/sL
@@ -0,0 +1 @@
+128 L 0
diff --git a/test/test_zsktrs_128_U/A b/test/test_zsktrs_128_U/A
new file mode 100644
index 0000000..077e405
Binary files /dev/null and b/test/test_zsktrs_128_U/A differ
diff --git a/test/test_zsktrs_128_U/B_orig b/test/test_zsktrs_128_U/B_orig
new file mode 100644
index 0000000..eb44548
Binary files /dev/null and b/test/test_zsktrs_128_U/B_orig differ
diff --git a/test/test_zsktrs_128_U/CMakeLists.txt b/test/test_zsktrs_128_U/CMakeLists.txt
new file mode 100644
index 0000000..284498a
--- /dev/null
+++ b/test/test_zsktrs_128_U/CMakeLists.txt
@@ -0,0 +1,8 @@
+get_filename_component(PARENT_DIR ${CMAKE_CURRENT_SOURCE_DIR} NAME)
+
+add_test_dependency_tree(
+  NAME "Test zsktrs wrapper U (128)"
+  COMMAND ${BASH_EXECUTABLE} cm.test.sh $<TARGET_FILE:test_zsktrs> U test_zsktrs.out
+  WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}
+  DEPENDENCY_TREE ${PARENT_DIR}
+  )
diff --git a/test/test_zsktrs_128_U/cm.test.sh b/test/test_zsktrs_128_U/cm.test.sh
new file mode 100755
index 0000000..f72b445
--- /dev/null
+++ b/test/test_zsktrs_128_U/cm.test.sh
@@ -0,0 +1,21 @@
+#!/bin/bash
+set -euo pipefail
+
+if [[ $# -gt 0 ]]; then
+	ZSKTRS=$1
+	UPLO=$2
+	OUT=$3
+else
+	source ../settings.sh
+fi
+
+if [ ! -f "$ZSKTRS" ]; then
+    echo "Executable $ZSKTRS does not exists"
+    exit 1
+fi
+
+echo " Test zsktrs wrapper " 
+cat s$UPLO | $ZSKTRS > $OUT
+[ $? -eq 0 ] && echo "Run without non-zero exit code" || exit 1
+
+grep OK $OUT
diff --git a/test/test_zsktrs_128_U/sU b/test/test_zsktrs_128_U/sU
new file mode 100644
index 0000000..95cd965
--- /dev/null
+++ b/test/test_zsktrs_128_U/sU
@@ -0,0 +1 @@
+128 U 0