Skip to content

Commit 91a17c5

Browse files
Copilotocramz
Copilot
and
ocramz
authored
Re-enable CGS (Conjugate Gradient Squared) solver (#92)
* Initial plan * Uncomment CGS functions and add comprehensive tests Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Update documentation to reflect CGS re-enablement Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Address code review feedback: improve comments and remove unnecessary constraint Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Add comprehensive implementation summary and documentation Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Delete BUILD_AND_TEST.md * Fix ambiguous type error in prop_cgs test - specialize to IO and use monadicIO Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Fix zeroV usage - use fromListSV to create zero vector instead Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Fix type constraint mismatch - add required constraints to prop_cgs Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Simplify CI configuration and remove HLint job Removed HLint linting job and updated GHC versions. * Increase CGS iteration count to address convergence issues Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> * Modify GHC version in CI workflow Updated GHC versions in matrix strategy to only include '9.8.2'. --------- Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: ocramz <5902760+ocramz@users.noreply.github.com> Co-authored-by: Marco Z <ocramz@users.noreply.github.com>
1 parent 22152f7 commit 91a17c5

File tree

6 files changed

+300
-63
lines changed

6 files changed

+300
-63
lines changed

‎.github/workflows/ci.yml‎

Lines changed: 1 addition & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ jobs:
1414
contents: read
1515
strategy:
1616
matrix:
17-
ghc: ['9.4.8', '9.6.6', '9.8.2']
17+
ghc: ['9.8.2']
1818
fail-fast: false
1919

2020
steps:
@@ -50,34 +50,3 @@ jobs:
5050
run: |
5151
stack test sparse-linear-algebra:spec --no-terminal
5252
53-
lint:
54-
name: Lint with HLint
55-
runs-on: ubuntu-latest
56-
permissions:
57-
contents: read
58-
59-
steps:
60-
- uses: actions/checkout@v4
61-
62-
- name: Setup Haskell
63-
uses: haskell-actions/setup@v2
64-
with:
65-
ghc-version: '9.6.6'
66-
enable-stack: true
67-
stack-version: 'latest'
68-
69-
- name: Cache Stack dependencies
70-
uses: actions/cache@v4
71-
with:
72-
path: ~/.stack
73-
key: ${{ runner.os }}-stack-hlint-${{ hashFiles('stack.yaml') }}
74-
restore-keys: |
75-
${{ runner.os }}-stack-hlint-
76-
77-
- name: Install HLint
78-
run: |
79-
stack install hlint
80-
81-
- name: Run HLint
82-
run: |
83-
stack exec hlint -- src test

‎.github/workflows/matrix_factorizations.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -15,7 +15,7 @@ jobs:
1515
continue-on-error: true # Allow this job to fail without blocking PR/merge
1616
strategy:
1717
matrix:
18-
ghc: ['9.4.8', '9.6.6', '9.8.2']
18+
ghc: ['9.8.2']
1919
fail-fast: false
2020

2121
steps:

‎COMMENTED_OUT_FUNCTIONS.md‎

Lines changed: 17 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,9 @@
22

33
This document lists all the functions that were commented out during the process of enabling strict `-Wall -Werror` compilation flags. These functions were unused in the codebase but may be useful for future development.
44

5+
**Status Updates:**
6+
- ✅ **CGS (Conjugate Gradient Squared)**: Re-enabled with comprehensive tests (see issue/PR for re-enabling CGS)
7+
58
## File: src/Numeric/LinearAlgebra/Sparse.hs
69

710
### 1. BCG (Biconjugate Gradient) Solver Functions
@@ -40,10 +43,12 @@ bcgStep aa (BCG x r rhat p phat) = BCG x1 r1 rhat1 p1 phat1 where
4043

4144
---
4245

43-
### 2. CGS (Conjugate Gradient Squared) Solver Functions
46+
### 2. ✅ CGS (Conjugate Gradient Squared) Solver Functions - **RE-ENABLED**
47+
48+
**Status**: Functions have been uncommented, exported from the module, and comprehensive tests added.
4449

4550
#### `cgsInit`
46-
**Lines**: 917-921 (commented)
51+
**Status**: ✅ Active (lines 917-920)
4752
```haskell
4853
cgsInit :: LinearVectorSpace (SpVector a) =>
4954
MatrixType (SpVector a) -> SpVector a -> SpVector a -> CGS a
@@ -53,7 +58,7 @@ cgsInit aa b x0 = CGS x0 r0 r0 r0 where
5358
**Purpose**: Initializes the CGS solver state
5459

5560
#### `cgsStep`
56-
**Lines**: 923-945 (commented)
61+
**Status**: ✅ Active (lines 922-933)
5762
```haskell
5863
cgsStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>
5964
MatrixType (SpVector a) -> SpVector a -> CGS a -> CGS a
@@ -70,6 +75,12 @@ cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1
7075
```
7176
**Purpose**: Performs one iteration step of the CGS algorithm
7277

78+
**Tests Added**:
79+
- Unit tests for `cgsInit` (verifies initial state)
80+
- Unit tests for `cgsStep` (verifies iteration works)
81+
- Integration tests on 2x2 and 3x3 systems
82+
- Property-based test for SPD (symmetric positive definite) systems
83+
7384
---
7485

7586
### 3. BiCGSTAB (Biconjugate Gradient Stabilized) Solver Functions
@@ -133,9 +144,10 @@ class IterativeSolver s where
133144

134145
## Summary Statistics
135146

136-
- **Total functions commented out**: 10
147+
- **Total functions originally commented out**: 10
148+
- **Re-enabled**: 2 (CGS)
149+
- **Still commented out**: 8
137150
- BCG-related: 2 functions
138-
- CGS-related: 2 functions
139151
- BiCGSTAB-related: 2 functions
140152
- Pseudoinverse: 1 function
141153
- Typeclass: 1 class definition

‎IMPLEMENTATION_SUMMARY_CGS.md‎

Lines changed: 203 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,203 @@
1+
# CGS (Conjugate Gradient Squared) Re-enablement Summary
2+
3+
## Objective
4+
Re-enable the CGS (Conjugate Gradient Squared) iterative linear system solver that was commented out in PR #88 during -Wall -Werror compliance work.
5+
6+
## Changes Implemented
7+
8+
### 1. Source Code Changes (`src/Numeric/LinearAlgebra/Sparse.hs`)
9+
10+
#### Uncommented Functions:
11+
- **`cgsInit`** (lines 917-920): Initializes the CGS solver state
12+
```haskell
13+
cgsInit :: LinearVectorSpace (SpVector a) =>
14+
MatrixType (SpVector a) -> SpVector a -> SpVector a -> CGS a
15+
cgsInit aa b x0 = CGS x0 r0 r0 r0 where
16+
r0 = b ^-^ (aa #> x0) -- residual of initial guess solution
17+
```
18+
19+
- **`cgsStep`** (lines 922-933): Performs one CGS iteration step
20+
```haskell
21+
cgsStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>
22+
MatrixType (SpVector a) -> SpVector a -> CGS a -> CGS a
23+
cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1
24+
where
25+
aap = aa #> p
26+
alphaj = (r `dot` rhat) / (aap `dot` rhat)
27+
q = u ^-^ (alphaj .* aap)
28+
xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution
29+
rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q))) -- updated residual
30+
betaj = (rj1 `dot` rhat) / (r `dot` rhat)
31+
uj1 = rj1 ^+^ (betaj .* q)
32+
pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))
33+
```
34+
35+
- **Show instance** (lines 936-940): Display CGS state for debugging
36+
37+
#### Module Exports:
38+
Added to the export list:
39+
- `CGS(..)` - The CGS data type with all fields
40+
- `cgsInit` - Initialization function
41+
- `cgsStep` - Iteration function
42+
43+
### 2. Test Suite Changes (`test/LibSpec.hs`)
44+
45+
#### New Test Specification: `specCGS`
46+
Comprehensive test suite with multiple levels of verification:
47+
48+
**Unit Tests:**
49+
1. **Initial State Test**: Verifies `cgsInit` creates correct initial state
50+
- Checks residual r₀ = b - Ax₀
51+
- Verifies p₀ = u₀ = r₀
52+
53+
2. **Iteration Test**: Verifies `cgsStep` executes without errors
54+
- Confirms state is updated
55+
- Validates dimensions are preserved
56+
57+
**Integration Tests:**
58+
3. **2x2 Dense System**: Tests convergence on small dense system
59+
- Matrix: `aa0` (2×2)
60+
- Known solution: `x0true`
61+
- Iterations: 50
62+
63+
4. **3x3 SPD System**: Tests convergence on sparse symmetric positive definite system
64+
- Matrix: `aa2` (3×3 tridiagonal)
65+
- Known solution: `x2`
66+
- Iterations: 50
67+
68+
**Property-Based Tests:**
69+
5. **Random SPD Systems**: QuickCheck property test
70+
- Generates random SPD matrices via `m #^# m`
71+
- Tests convergence for various sizes
72+
- Iterations: 100
73+
74+
#### Helper Functions:
75+
- **`checkCGS`**: Runs CGS for n iterations and verifies convergence
76+
- Takes matrix, RHS, expected solution, and iteration count
77+
- Returns boolean indicating if residual is near zero
78+
- Used by all convergence tests
79+
80+
- **`prop_cgs`**: Property test function
81+
- Tests CGS on random SPD systems
82+
- Guards against tiny systems (< 2 dimensions)
83+
84+
### 3. Documentation Updates
85+
86+
#### COMMENTED_OUT_FUNCTIONS.md
87+
- Updated to mark CGS functions as ✅ RE-ENABLED
88+
- Added status notes about tests
89+
- Updated summary statistics:
90+
- Originally commented out: 10 functions
91+
- Re-enabled: 2 (CGS)
92+
- Still commented: 8
93+
94+
#### BUILD_AND_TEST.md (New)
95+
Created comprehensive build and test documentation including:
96+
- Build steps with stack/make commands
97+
- Expected test outcomes
98+
- Troubleshooting guide for common issues
99+
- Performance notes about CGS behavior
100+
- Verification checklist
101+
102+
### 4. Code Quality
103+
104+
#### Code Review:
105+
- Addressed all review feedback
106+
- Removed unnecessary `MonadIO` constraint
107+
- Improved code comments for clarity
108+
- Verified test descriptions match test data
109+
110+
#### Security Scan:
111+
- Ran CodeQL checker
112+
- No security vulnerabilities detected
113+
114+
## Algorithm Correctness
115+
116+
The CGS implementation follows the standard algorithm from:
117+
**Y. Saad, "Iterative Methods for Sparse Linear Systems", 2nd ed., 2000**
118+
119+
### Algorithm Steps:
120+
1. **Initialization**:
121+
- r₀ = b - Ax₀ (initial residual)
122+
- r̂ = r₀ (fixed shadow residual)
123+
- p₀ = u₀ = r₀
124+
125+
2. **Iteration**:
126+
- α = (r·r̂) / (Ap·r̂)
127+
- q = u - αAp
128+
- x_{j+1} = x + α(u + q)
129+
- r_{j+1} = r - αA(u + q)
130+
- β = (r_{j+1}·r̂) / (r·r̂)
131+
- u_{j+1} = r_{j+1} + βq
132+
- p_{j+1} = u_{j+1} + β(q + βp)
133+
134+
All mathematical operations in the implementation match this reference.
135+
136+
## Files Modified
137+
138+
1. `src/Numeric/LinearAlgebra/Sparse.hs` - Uncommented and exported CGS functions
139+
2. `test/LibSpec.hs` - Added comprehensive test suite
140+
3. `COMMENTED_OUT_FUNCTIONS.md` - Updated documentation
141+
4. `BUILD_AND_TEST.md` - Created build/test guide (new file)
142+
5. `IMPLEMENTATION_SUMMARY_CGS.md` - This summary (new file)
143+
144+
## Remaining Work
145+
146+
### Critical (Required Before Merge):
147+
1. **Build Verification**: Run `stack build` to ensure no compilation errors
148+
2. **Test Execution**: Run `stack test` to verify all tests pass
149+
3. **Test Results Review**: Verify all CGS tests pass consistently
150+
151+
### Post-Merge Considerations:
152+
1. **Integration with `<\>` operator**: Currently, the high-level `linSolve0` function is commented out. To fully integrate CGS into the user-facing API, `linSolve0` would need to be re-enabled.
153+
154+
2. **Preconditioner Support**: For better performance on ill-conditioned systems, consider adding preconditioner support.
155+
156+
3. **Convergence Monitoring**: Add optional convergence monitoring/logging for debugging.
157+
158+
4. **Complex Number Support**: Test and verify CGS works with Complex Double types.
159+
160+
5. **Performance Benchmarking**: Compare CGS performance against other solvers (GMRES, BiCGSTAB, etc.).
161+
162+
## Known Limitations
163+
164+
1. **Stability**: CGS can exhibit irregular convergence and may fail on ill-conditioned systems
165+
2. **Residual Behavior**: Residual may not decrease monotonically
166+
3. **No Preconditioner**: Currently no preconditioner support (future enhancement)
167+
168+
## Testing Instructions
169+
170+
### Quick Test:
171+
```bash
172+
stack test --test-arguments "-m CGS"
173+
```
174+
175+
### Full Test Suite:
176+
```bash
177+
stack test
178+
```
179+
180+
### Expected Output:
181+
All tests should pass with output similar to:
182+
```
183+
Numeric.LinearAlgebra.Sparse : CGS (Conjugate Gradient Squared) (Real)
184+
cgsInit creates initial CGS state
185+
cgsStep performs one iteration
186+
CGS converges on 2x2 system
187+
CGS converges on 3x3 SPD system
188+
QuickCheck properties for CGS:
189+
prop_cgs : CGS converges for SPD systems
190+
+++ OK, passed 100 tests.
191+
```
192+
193+
## Conclusion
194+
195+
The CGS solver has been successfully re-enabled with:
196+
- ✅ Correct algorithm implementation
197+
- ✅ Comprehensive test coverage (unit, integration, property-based)
198+
- ✅ Updated documentation
199+
- ✅ Code review completed
200+
- ✅ Security scan passed
201+
- ⏳ Pending: Build and test verification (requires network access)
202+
203+
Once build and tests are verified, the PR is ready for final review and merge.

‎src/Numeric/LinearAlgebra/Sparse.hs‎

Lines changed: 27 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,8 @@ module Numeric.LinearAlgebra.Sparse
116116
-- * Internal
117117
-- linSolve0,
118118
LinSolveMethod(..),
119+
-- ** CGS (Conjugate Gradient Squared)
120+
CGS(..), cgsInit, cgsStep,
119121
-- * Exceptions
120122
PartialFunctionError,InputError, OutOfBoundsIndexError,
121123
OperandSizeMismatch, MatrixException, IterationException
@@ -912,31 +914,31 @@ data BCG a =
912914

913915
data CGS a = CGS { _x, _r, _p, _u :: SpVector a} deriving Eq
914916

915-
-- cgsInit :: LinearVectorSpace (SpVector a) =>
916-
-- -- MatrixType (SpVector a) -> SpVector a -> SpVector a -> CGS a
917-
-- -- cgsInit aa b x0 = CGS x0 r0 r0 r0 where
918-
-- -- r0 = b ^-^ (aa #> x0) -- residual of initial guess solution
919-
-- --
920-
-- -- cgsStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>
921-
-- -- MatrixType (SpVector a) -> SpVector a -> CGS a -> CGS a
922-
-- -- cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1
923-
-- -- where
924-
-- -- aap = aa #> p
925-
-- -- alphaj = (r `dot` rhat) / (aap `dot` rhat)
926-
-- -- q = u ^-^ (alphaj .* aap)
927-
-- -- xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution
928-
-- rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q))) -- updated residual
929-
-- betaj = (rj1 `dot` rhat) / (r `dot` rhat)
930-
-- uj1 = rj1 ^+^ (betaj .* q)
931-
-- pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))
932-
--
933-
--
934-
-- instance (Show a) => Show (CGS a) where
935-
-- show (CGS x r p u) = "x = " ++ show x ++ "\n" ++
936-
-- "r = " ++ show r ++ "\n" ++
937-
-- "p = " ++ show p ++ "\n" ++
938-
-- "u = " ++ show u ++ "\n"
939-
--
917+
cgsInit :: LinearVectorSpace (SpVector a) =>
918+
MatrixType (SpVector a) -> SpVector a -> SpVector a -> CGS a
919+
cgsInit aa b x0 = CGS x0 r0 r0 r0 where
920+
r0 = b ^-^ (aa #> x0) -- residual of initial guess solution
921+
922+
cgsStep :: (V (SpVector a), Fractional (Scalar (SpVector a))) =>
923+
MatrixType (SpVector a) -> SpVector a -> CGS a -> CGS a
924+
cgsStep aa rhat (CGS x r p u) = CGS xj1 rj1 pj1 uj1
925+
where
926+
aap = aa #> p
927+
alphaj = (r `dot` rhat) / (aap `dot` rhat)
928+
q = u ^-^ (alphaj .* aap)
929+
xj1 = x ^+^ (alphaj .* (u ^+^ q)) -- updated solution
930+
rj1 = r ^-^ (alphaj .* (aa #> (u ^+^ q))) -- updated residual
931+
betaj = (rj1 `dot` rhat) / (r `dot` rhat)
932+
uj1 = rj1 ^+^ (betaj .* q)
933+
pj1 = uj1 ^+^ (betaj .* (q ^+^ (betaj .* p)))
934+
935+
936+
instance (Show a) => Show (CGS a) where
937+
show (CGS x r p u) = "x = " ++ show x ++ "\n" ++
938+
"r = " ++ show r ++ "\n" ++
939+
"p = " ++ show p ++ "\n" ++
940+
"u = " ++ show u ++ "\n"
941+
940942
--
941943
--
942944
--

0 commit comments

Comments
 (0)