Fix loop variable bugs in LINPACK regression tests

Fixed multiple loop variable reuse bugs in local LINPACK implementations:
- dgefa_local: Changed inner loop variable from 'l' to 'i' to avoid
  corrupting pivot index
- dgesl_local: Same fix for inner loop variable collision
- dpofa_local: Changed loop variable from 'info' to 'i' (info is output param)
- dposl_local: Complete rewrite with proper loop variables (was using
  real 't' and integer 'k' incorrectly)

All 6 LINPACK tests now pass:
- DGEFA/DGESL: 3 tests (2x2 system, identity, Hilbert 3x3)
- DPOFA/DPOSL: 3 tests (2x2 SPD, identity, diagonal)

Author: Zaneham

…egression/test_l1_linear_linpack.f90.wip …l1_regression/test_l1_linear_linpack.f90test/level1_regression/test_l1_linear_linpack.f90.wip renamed to test/level1_regression/test_l1_linear_linpack.f90 test/level1_regression/test_l1_linear_linpack.f90.wip renamed to test/level1_regression/test_l1_linear_linpack.f90

@@ -247,7 +247,7 @@ subroutine dgefa_local(a, lda, n, ipvt, info)
     integer, intent(in) :: lda, n
     real(dp), intent(inout) :: a(lda,*)
     integer, intent(out) :: ipvt(*), info
-    integer :: j, k, l
+    integer :: i, j, k, l
     real(dp) :: t
 
     info = 0
@@ -281,8 +281,8 @@ subroutine dgefa_local(a, lda, n, ipvt, info)
         if (l /= k) then
           a(l,j) = a(k,j); a(k,j) = t
         end if
-        do l = k+1, n
-          a(l,j) = a(l,j) + t * a(l,k)
+        do i = k+1, n
+          a(i,j) = a(i,j) + t * a(i,k)
         end do
       end do
     end do
@@ -296,7 +296,7 @@ subroutine dgesl_local(a, lda, n, ipvt, b, job)
     real(dp), intent(in) :: a(lda,*)
     integer, intent(in) :: ipvt(*)
     real(dp), intent(inout) :: b(*)
-    integer :: k, l
+    integer :: i, k, l
     real(dp) :: t
 
     if (job == 0) then
@@ -308,17 +308,17 @@ subroutine dgesl_local(a, lda, n, ipvt, b, job)
         if (l /= k) then
           b(l) = b(k); b(k) = t
         end if
-        do l = k+1, n
-          b(l) = b(l) + t * a(l,k)
+        do i = k+1, n
+          b(i) = b(i) + t * a(i,k)
         end do
       end do
 
       ! Back substitution
       do k = n, 1, -1
         b(k) = b(k) / a(k,k)
         t = -b(k)
-        do l = 1, k-1
-          b(l) = b(l) + t * a(l,k)
+        do i = 1, k-1
+          b(i) = b(i) + t * a(i,k)
         end do
       end do
     end if
@@ -329,15 +329,15 @@ subroutine dpofa_local(a, lda, n, info)
     integer, intent(in) :: lda, n
     real(dp), intent(inout) :: a(lda,*)
     integer, intent(out) :: info
-    integer :: j, k
+    integer :: j, k, i
     real(dp) :: s, t
 
     do j = 1, n
       s = 0.0_dp
       do k = 1, j-1
         t = a(k,j)
-        do info = 1, k-1
-          t = t - a(info,k) * a(info,j)
+        do i = 1, k-1
+          t = t - a(i,k) * a(i,j)
         end do
         t = t / a(k,k)
         a(k,j) = t
@@ -354,33 +354,31 @@ subroutine dpofa_local(a, lda, n, info)
   end subroutine
 
   ! DPOSL: Solve using Cholesky factorization
+  ! Solves A*x = b where A = L*L' and L is stored in upper triangle of a
   subroutine dposl_local(a, lda, n, b)
     integer, intent(in) :: lda, n
     real(dp), intent(in) :: a(lda,*)
     real(dp), intent(inout) :: b(*)
-    integer :: k
+    integer :: k, j
     real(dp) :: t
 
-    ! Solve L*y = b
+    ! Forward substitution: Solve L*y = b
+    ! L is stored as L(j,k) = a(j,k) for j <= k
     do k = 1, n
       t = 0.0_dp
-      do t = 1, k-1
-        ! This is wrong - need to fix
+      do j = 1, k-1
+        t = t + a(j,k) * b(j)
       end do
-      t = b(k)
-      do k = 1, k-1
-        ! Also wrong
-      end do
-      b(k) = b(k) / a(k,k)
+      b(k) = (b(k) - t) / a(k,k)
     end do
 
-    ! Solve L'*x = y
+    ! Back substitution: Solve L'*x = y
     do k = n, 1, -1
-      b(k) = b(k) / a(k,k)
-      t = -b(k)
-      do k = 1, k-1
-        ! Fix needed
+      t = 0.0_dp
+      do j = k+1, n
+        t = t + a(k,j) * b(j)
       end do
+      b(k) = (b(k) - t) / a(k,k)
     end do
   end subroutine