Skip to content

Commit fd200ce

Browse files
sriharikrishnasriharikrishna
committed
Outlining of banded matrices and vector codes
1 parent 72655f5 commit fd200ce

312 files changed

Lines changed: 18074 additions & 15610 deletions

File tree

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎BLAS/Makefile‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -297,31 +297,31 @@ $(BUILD_DIR)/libdiffblas_d.a: compile-d $(DIFFSIZES_ACCESS_OBJ)
297297
@echo "Created libdiffblas_d.a with $$(ls $(BUILD_DIR)/*_d.o 2>/dev/null | wc -w) objects"
298298

299299
$(BUILD_DIR)/libdiffblas_d.so: compile-d
300-
@$(FC) -shared -o $@ $$(ls $(BUILD_DIR)/*_d.o 2>/dev/null)
300+
@objs="$$(ls $(BUILD_DIR)/*_d.o 2>/dev/null)"; if [ -n "$$objs" ]; then $(FC) -shared -o $@ $$objs; else touch $@; fi
301301

302302
# Single library for all reverse mode differentiated code
303303
$(BUILD_DIR)/libdiffblas_b.a: compile-b $(DIFFSIZES_ACCESS_OBJ)
304304
@ar rcs $@ $$(ls $(BUILD_DIR)/*_b.o 2>/dev/null) $(BUILD_DIR)/adStack.o $(DIFFSIZES_ACCESS_OBJ)
305305
@echo "Created libdiffblas_b.a with $$(ls $(BUILD_DIR)/*_b.o 2>/dev/null | wc -w) objects"
306306

307307
$(BUILD_DIR)/libdiffblas_b.so: compile-b $(DIFFSIZES_ACCESS_OBJ)
308-
@$(FC) -shared -o $@ $$(ls $(BUILD_DIR)/*_b.o 2>/dev/null) $(BUILD_DIR)/adStack.o $(DIFFSIZES_ACCESS_OBJ)
308+
@objs="$$(ls $(BUILD_DIR)/*_b.o 2>/dev/null)"; if [ -n "$$objs" ]; then $(FC) -shared -o $@ $$objs $(BUILD_DIR)/adStack.o $(DIFFSIZES_ACCESS_OBJ); else touch $@; fi
309309

310310
# Single library for all vector forward mode differentiated code
311311
$(BUILD_DIR)/libdiffblas_dv.a: compile-dv $(DIFFSIZES_ACCESS_OBJ)
312312
@ar rcs $@ $$(ls $(BUILD_DIR)/*_dv.o 2>/dev/null) $(BUILD_DIR)/DIFFSIZES.o $(DIFFSIZES_ACCESS_OBJ)
313313
@echo "Created libdiffblas_dv.a with $$(ls $(BUILD_DIR)/*_dv.o 2>/dev/null | wc -w) objects"
314314

315315
$(BUILD_DIR)/libdiffblas_dv.so: compile-dv
316-
@$(FC) -shared -o $@ $$(ls $(BUILD_DIR)/*_dv.o 2>/dev/null) $(BUILD_DIR)/DIFFSIZES.o
316+
@objs="$$(ls $(BUILD_DIR)/*_dv.o 2>/dev/null)"; if [ -n "$$objs" ]; then $(FC) -shared -o $@ $$objs $(BUILD_DIR)/DIFFSIZES.o; else touch $@; fi
317317

318318
# Single library for all vector reverse mode differentiated code
319319
$(BUILD_DIR)/libdiffblas_bv.a: compile-bv $(DIFFSIZES_ACCESS_OBJ)
320320
@ar rcs $@ $$(ls $(BUILD_DIR)/*_bv.o 2>/dev/null) $(BUILD_DIR)/adStack.o $(BUILD_DIR)/DIFFSIZES.o $(DIFFSIZES_ACCESS_OBJ)
321321
@echo "Created libdiffblas_bv.a with $$(ls $(BUILD_DIR)/*_bv.o 2>/dev/null | wc -w) objects"
322322

323323
$(BUILD_DIR)/libdiffblas_bv.so: compile-bv $(DIFFSIZES_ACCESS_OBJ)
324-
@$(FC) -shared -o $@ $$(ls $(BUILD_DIR)/*_bv.o 2>/dev/null) $(BUILD_DIR)/adStack.o $(BUILD_DIR)/DIFFSIZES.o $(DIFFSIZES_ACCESS_OBJ)
324+
@objs="$$(ls $(BUILD_DIR)/*_bv.o 2>/dev/null)"; if [ -n "$$objs" ]; then $(FC) -shared -o $@ $$objs $(BUILD_DIR)/adStack.o $(BUILD_DIR)/DIFFSIZES.o $(DIFFSIZES_ACCESS_OBJ); else touch $@; fi
325325

326326
# Note: Original BLAS functions come from $(BLAS_LIB) (librefblas in LAPACKDIR)
327327
# No need to build a separate liborigblas

‎BLAS/test/test_caxpy.f90‎

Lines changed: 6 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -108,17 +108,17 @@ subroutine run_test_for_size(n, passed)
108108
write(*,*) 'Function calls completed successfully'
109109

110110
! Numerical differentiation check
111-
call check_derivatives_numerically(n, nsize, cx_orig, cy_orig, ca_orig, cx_d_orig, cy_d_orig, ca_d_orig, cy_d, passed)
111+
call check_derivatives_numerically(n, nsize, ca_orig, cx_orig, cy_orig, ca_d_orig, cx_d_orig, cy_d_orig, cy_d, passed)
112112

113113
end subroutine run_test_for_size
114114

115-
subroutine check_derivatives_numerically(n, nsize, cx_orig, cy_orig, ca_orig, cx_d_orig, cy_d_orig, ca_d_orig, cy_d, passed)
115+
subroutine check_derivatives_numerically(n, nsize, ca_orig, cx_orig, cy_orig, ca_d_orig, cx_d_orig, cy_d_orig, cy_d, passed)
116116
implicit none
117117
integer, intent(in) :: n
118118
integer, intent(in) :: nsize
119+
complex(4), intent(in) :: ca_orig, ca_d_orig
119120
complex(4), intent(in) :: cx_orig(n), cx_d_orig(n)
120121
complex(4), intent(in) :: cy_orig(n), cy_d_orig(n)
121-
complex(4), intent(in) :: ca_orig, ca_d_orig
122122
complex(4), intent(in) :: cy_d(n)
123123
logical, intent(out) :: passed
124124

@@ -129,9 +129,9 @@ subroutine check_derivatives_numerically(n, nsize, cx_orig, cy_orig, ca_orig, cx
129129
logical :: has_large_errors
130130
complex(4), dimension(n) :: cy_forward, cy_backward
131131
integer :: i, j
132+
complex(4) :: ca
132133
complex(4), dimension(n) :: cx
133134
complex(4), dimension(n) :: cy
134-
complex(4) :: ca
135135

136136
max_error = 0.0e0
137137
has_large_errors = .false.
@@ -140,16 +140,16 @@ subroutine check_derivatives_numerically(n, nsize, cx_orig, cy_orig, ca_orig, cx
140140
write(*,*) 'Step size h =', h
141141

142142
! Forward perturbation: f(x + h)
143+
ca = ca_orig + h * ca_d_orig
143144
cx = cx_orig + h * cx_d_orig
144145
cy = cy_orig + h * cy_d_orig
145-
ca = ca_orig + h * ca_d_orig
146146
call caxpy(nsize, ca, cx, 1, cy, 1)
147147
cy_forward = cy
148148

149149
! Backward perturbation: f(x - h)
150+
ca = ca_orig - h * ca_d_orig
150151
cx = cx_orig - h * cx_d_orig
151152
cy = cy_orig - h * cy_d_orig
152-
ca = ca_orig - h * ca_d_orig
153153
call caxpy(nsize, ca, cx, 1, cy, 1)
154154
cy_backward = cy
155155

‎BLAS/test/test_caxpy_vector_forward.f90‎

Lines changed: 67 additions & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -46,78 +46,86 @@ program test_caxpy_vector_forward
4646
n = test_sizes(itest)
4747
write(*,*) 'Testing CAXPY (Vector Forward, n =', n, ')'
4848

49-
! Initialize test parameters
50-
nsize = n
51-
incx_val = 1
52-
incy_val = 1
49+
call run_test_for_size(n, passed)
50+
all_passed = all_passed .and. passed
51+
end do
52+
if (all_passed) then
53+
write(*,*) 'PASS: Vector forward mode - all sizes completed successfully'
54+
else
55+
write(*,*) 'FAIL: Vector forward mode - one or more sizes had derivative errors'
56+
end if
5357

54-
! Initialize test data with random numbers
55-
! Initialize random seed for reproducible results
56-
seed_array = 42
57-
call random_seed(put=seed_array)
58+
contains
5859

59-
call random_number(temp_real)
60-
call random_number(temp_imag)
61-
ca = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
62-
do i = 1, max_size
63-
call random_number(temp_real)
64-
call random_number(temp_imag)
65-
cx(i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
66-
end do
67-
do i = 1, max_size
68-
call random_number(temp_real)
69-
call random_number(temp_imag)
70-
cy(i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
71-
end do
60+
subroutine run_test_for_size(n, passed)
61+
implicit none
62+
integer, intent(in) :: n
63+
logical, intent(out) :: passed
7264

73-
! Initialize input derivatives to random values (exactly like scalar mode)
74-
do idir = 1, nbdirs
65+
! Initialize test parameters
66+
nsize = n
67+
incx_val = 1
68+
incy_val = 1
69+
70+
! Initialize test data with random numbers
71+
! Initialize random seed for reproducible results
72+
seed_array = 42
73+
call random_seed(put=seed_array)
74+
7575
call random_number(temp_real)
7676
call random_number(temp_imag)
77-
ca_dv(idir) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
78-
end do
79-
do idir = 1, nbdirs
77+
ca = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
8078
do i = 1, max_size
8179
call random_number(temp_real)
8280
call random_number(temp_imag)
83-
cx_dv(idir,i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
81+
cx(i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
8482
end do
85-
end do
86-
do idir = 1, nbdirs
8783
do i = 1, max_size
8884
call random_number(temp_real)
8985
call random_number(temp_imag)
90-
cy_dv(idir,i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
86+
cy(i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
9187
end do
92-
end do
93-
94-
write(*,*) 'Testing CAXPY (Vector Forward Mode)'
95-
! Store original values before any function calls (critical for INOUT parameters)
96-
ca_orig = ca
97-
ca_dv_orig = ca_dv
98-
cx_orig = cx
99-
cx_dv_orig = cx_dv
100-
cy_orig = cy
101-
cy_dv_orig = cy_dv
102-
103-
! Call the vector mode differentiated function
104-
105-
call caxpy_dv(nsize, ca, ca_dv, cx, cx_dv, incx_val, cy, cy_dv, incy_val, nbdirs)
106-
107-
! Print results and compare
108-
write(*,*) 'Function calls completed successfully'
109-
110-
! Numerical differentiation check
111-
call check_derivatives_numerically(passed)
112-
all_passed = all_passed .and. passed
113-
end do
114-
if (all_passed) then
115-
write(*,*) 'PASS: Vector forward mode - all sizes completed successfully'
116-
else
117-
write(*,*) 'FAIL: Vector forward mode - one or more sizes had derivative errors'
118-
end if
119-
120-
contains
88+
89+
! Initialize input derivatives to random values (exactly like scalar mode)
90+
do idir = 1, nbdirs
91+
call random_number(temp_real)
92+
call random_number(temp_imag)
93+
ca_dv(idir) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
94+
end do
95+
do idir = 1, nbdirs
96+
do i = 1, max_size
97+
call random_number(temp_real)
98+
call random_number(temp_imag)
99+
cx_dv(idir,i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
100+
end do
101+
end do
102+
do idir = 1, nbdirs
103+
do i = 1, max_size
104+
call random_number(temp_real)
105+
call random_number(temp_imag)
106+
cy_dv(idir,i) = cmplx(temp_real, temp_imag) * (2.0,2.0) - (1.0,1.0)
107+
end do
108+
end do
109+
110+
write(*,*) 'Testing CAXPY (Vector Forward Mode)'
111+
! Store original values before any function calls (critical for INOUT parameters)
112+
ca_orig = ca
113+
ca_dv_orig = ca_dv
114+
cx_orig = cx
115+
cx_dv_orig = cx_dv
116+
cy_orig = cy
117+
cy_dv_orig = cy_dv
118+
119+
! Call the vector mode differentiated function
120+
121+
call caxpy_dv(nsize, ca, ca_dv, cx, cx_dv, incx_val, cy, cy_dv, incy_val, nbdirs)
122+
123+
! Print results and compare
124+
write(*,*) 'Function calls completed successfully'
125+
126+
! Numerical differentiation check
127+
call check_derivatives_numerically(passed)
128+
end subroutine run_test_for_size
121129

122130
subroutine check_derivatives_numerically(passed)
123131
implicit none

‎BLAS/test/test_caxpy_vector_reverse.f90‎

Lines changed: 62 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -59,59 +59,7 @@ program test_caxpy_vector_reverse
5959
n = test_sizes(itest)
6060
write(*,*) 'Testing CAXPY (Vector Reverse, n =', n, ')'
6161

62-
! Initialize primal values
63-
nsize = n
64-
call random_number(temp_real)
65-
call random_number(temp_imag)
66-
ca = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
67-
do i = 1, n
68-
call random_number(temp_real)
69-
call random_number(temp_imag)
70-
cx(i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
71-
end do
72-
incx_val = 1
73-
do i = 1, n
74-
call random_number(temp_real)
75-
call random_number(temp_imag)
76-
cy(i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
77-
end do
78-
incy_val = 1
79-
80-
! Store original primal values
81-
ca_orig = ca
82-
cx_orig = cx
83-
cy_orig = cy
84-
85-
! Initialize output adjoints (cotangents) with random values for each direction
86-
! These are the 'seeds' for reverse mode
87-
do k = 1, nbdirs
88-
do i = 1, n
89-
call random_number(temp_real)
90-
call random_number(temp_imag)
91-
cyb(k,i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
92-
end do
93-
end do
94-
95-
! Initialize input adjoints to zero (they will be computed)
96-
! Note: Inout parameters are skipped - they already have output adjoints initialized
97-
cab = 0.0
98-
cxb = 0.0
99-
100-
! Save original cotangent seeds for OUTPUT/INOUT parameters (before function call)
101-
cyb_orig = cyb
102-
103-
! Set ISIZE globals required by differentiated routine (dimension 2 of arrays).
104-
! ISIZE1OF* (vectors): use n to match adjoint array size; ISIZE2OF* (matrices): use max_size.
105-
call set_ISIZE1OFCx(n)
106-
107-
! Call reverse vector mode differentiated function
108-
call caxpy_bv(nsize, ca, cab, cx, cxb, incx_val, cy, cyb, incy_val, nbdirs)
109-
110-
! Reset ISIZE globals to uninitialized (-1) for completeness
111-
call set_ISIZE1OFCx(-1)
112-
113-
! VJP Verification using finite differences
114-
call check_vjp_numerically(passed)
62+
call run_test_for_size(n, passed)
11563
all_passed = all_passed .and. passed
11664
end do
11765
if (all_passed) then
@@ -122,6 +70,66 @@ program test_caxpy_vector_reverse
12270

12371
contains
12472

73+
subroutine run_test_for_size(n, passed)
74+
implicit none
75+
integer, intent(in) :: n
76+
logical, intent(out) :: passed
77+
78+
! Initialize primal values
79+
nsize = n
80+
call random_number(temp_real)
81+
call random_number(temp_imag)
82+
ca = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
83+
do i = 1, n
84+
call random_number(temp_real)
85+
call random_number(temp_imag)
86+
cx(i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
87+
end do
88+
incx_val = 1
89+
do i = 1, n
90+
call random_number(temp_real)
91+
call random_number(temp_imag)
92+
cy(i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
93+
end do
94+
incy_val = 1
95+
96+
! Store original primal values
97+
ca_orig = ca
98+
cx_orig = cx
99+
cy_orig = cy
100+
101+
! Initialize output adjoints (cotangents) with random values for each direction
102+
! These are the 'seeds' for reverse mode
103+
do k = 1, nbdirs
104+
do i = 1, n
105+
call random_number(temp_real)
106+
call random_number(temp_imag)
107+
cyb(k,i) = cmplx(temp_real * 2.0 - 1.0, temp_imag * 2.0 - 1.0)
108+
end do
109+
end do
110+
111+
! Initialize input adjoints to zero (they will be computed)
112+
! Note: Inout parameters are skipped - they already have output adjoints initialized
113+
cab = 0.0
114+
cxb = 0.0
115+
116+
! Save original cotangent seeds for OUTPUT/INOUT parameters (before function call)
117+
cyb_orig = cyb
118+
119+
! Set ISIZE globals required by differentiated routine (dimension 2 of arrays).
120+
! ISIZE1OF* (vectors): use n to match adjoint array size; ISIZE2OF* (matrices): use max_size.
121+
call set_ISIZE1OFCx(n)
122+
123+
! Call reverse vector mode differentiated function
124+
call caxpy_bv(nsize, ca, cab, cx, cxb, incx_val, cy, cyb, incy_val, nbdirs)
125+
126+
! Reset ISIZE globals to uninitialized (-1) for completeness
127+
call set_ISIZE1OFCx(-1)
128+
129+
! VJP Verification using finite differences
130+
call check_vjp_numerically(passed)
131+
end subroutine run_test_for_size
132+
125133
subroutine check_vjp_numerically(passed)
126134
implicit none
127135
logical, intent(out) :: passed
@@ -196,6 +204,7 @@ subroutine check_vjp_numerically(passed)
196204
! For INOUT parameters: use cb directly (it contains the computed input adjoint after reverse pass)
197205
! For pure inputs: use adjoint directly
198206
vjp_ad = 0.0
207+
vjp_ad = vjp_ad + real(conjg(ca_dir) * cab(k))
199208
! Compute and sort products for cx
200209
n_products = n
201210
do i = 1, n
@@ -214,7 +223,6 @@ subroutine check_vjp_numerically(passed)
214223
do i = 1, n_products
215224
vjp_ad = vjp_ad + temp_products(i)
216225
end do
217-
vjp_ad = vjp_ad + real(conjg(ca_dir) * cab(k))
218226

219227
! Error check: |vjp_fd - vjp_ad| > atol + rtol * |vjp_ad|
220228
abs_error = abs(vjp_fd - vjp_ad)

0 commit comments

Comments
 (0)