Skip to content

Commit bd80504

Browse files
committed
Added test for a tall matrix with rank deficiency.
1 parent dd9b022 commit bd80504

File tree

1 file changed

+81
-0
lines changed

1 file changed

+81
-0
lines changed

test/linalg/test_linalg_pivoting_qr.fypp

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,87 @@ module test_linalg_pivoting_qr
112112
if (allocated(error)) return
113113
end block
114114

115+
!----------------------------------------------------
116+
!----- Tall matrix with rank deficiency -----
117+
!----------------------------------------------------
118+
block
119+
integer(ilp), parameter :: m = 15_ilp
120+
integer(ilp), parameter :: n = 4_ilp
121+
integer(ilp), parameter :: k = min(m,n)
122+
real(${rk}$), parameter :: tol = 100*sqrt(epsilon(0.0_${rk}$))
123+
${rt}$ :: a(m,n),aorig(m,n),q(m,m),r(m,n),qred(m,k),rred(k,n),qerr(m,6),rerr(6,n)
124+
real(${rk}$) :: rea(m,n),ima(m,n)
125+
integer(ilp) :: pivots(n), i, j
126+
integer(ilp) :: lwork
127+
${rt}$, allocatable :: work(:)
128+
type(linalg_state_type) :: state
129+
130+
call random_number(rea)
131+
#:if rt.startswith('complex')
132+
call random_number(ima)
133+
a = cmplx(rea,ima,kind=${rk}$)
134+
#:else
135+
a = rea
136+
#:endif
137+
a(:, 3) = 0.0_${rk}$ ! Zero-out column to test rank-deficiency.
138+
aorig = a
139+
140+
! 1) QR factorization with full matrices. Input NaNs to be sure Q and R are OK on return
141+
q = ieee_value(0.0_${rk}$,ieee_quiet_nan)
142+
r = ieee_value(0.0_${rk}$,ieee_quiet_nan)
143+
call qr(a, q, r, pivots, err=state)
144+
145+
! Check return code
146+
call check(error,state%ok(),state%print())
147+
if (allocated(error)) return
148+
149+
! Check solution
150+
call check(error, all(abs(a(:, pivots)-matmul(q,r))<tol), 'converged solution (fulle)')
151+
if (allocated(error)) return
152+
153+
! 2) QR factorization with reduced matrices
154+
call qr(a, qred, rred, pivots, err=state)
155+
156+
! Check return code
157+
call check(error,state%ok(),state%print())
158+
if (allocated(error)) return
159+
160+
! Check solution
161+
call check(error, all(abs(a(:, pivots)-matmul(qred,rred))<tol), 'converged solution (reduced)')
162+
if (allocated(error)) return
163+
164+
! 3) overwrite A
165+
call qr(a, qred, rred, pivots, overwrite_a=.true., err=state)
166+
167+
! Check return code
168+
call check(error,state%ok(),state%print())
169+
if (allocated(error)) return
170+
171+
! Check solution
172+
call check(error, all(abs(aorig(:, pivots)-matmul(qred,rred))<tol), 'converged solution (overwrite A)')
173+
if (allocated(error)) return
174+
175+
! 4) External storage option
176+
a = aorig
177+
call qr_space(a, lwork, pivoting=.true.)
178+
allocate(work(lwork))
179+
call qr(a, q, r, pivots, storage=work, err=state)
180+
181+
! Check return code
182+
call check(error,state%ok(),state%print())
183+
if (allocated(error)) return
184+
185+
! Check solution
186+
call check(error, all(abs(a(:, pivots)-matmul(q,r))<tol), 'converged solution (external storage)')
187+
if (allocated(error)) return
188+
189+
! Check that an invalid problem size returns an error
190+
a = aorig
191+
call qr(a, qerr, rerr, pivots, err=state)
192+
call check(error,state%error(),'invalid matrix sizes')
193+
if (allocated(error)) return
194+
end block
195+
115196
!-------------------------------
116197
!----- Wide matrix -----
117198
!-------------------------------

0 commit comments

Comments
 (0)