Skip to content

Commit 6c7e38b

Browse files
thijssteelthijssteel
committed
add rotc to blas
1 parent 6061809 commit 6c7e38b

File tree

7 files changed

+1534
-4
lines changed

7 files changed

+1534
-4
lines changed

‎BLAS/SRC/CMakeLists.txt‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -82,15 +82,15 @@ set(ZBLAS2 zgemv.f zgbmv.f zhemv.f zhbmv.f zhpmv.f
8282
#---------------------------------------------------------
8383
# Level 3 BLAS
8484
#---------------------------------------------------------
85-
set(SBLAS3 sgemm.f ssymm.f ssyrk.f ssyr2k.f strmm.f strsm.f sgemmtr.f)
85+
set(SBLAS3 sgemm.f ssymm.f ssyrk.f ssyr2k.f strmm.f strsm.f sgemmtr.f srotc.f90)
8686

8787
set(CBLAS3 cgemm.f csymm.f csyrk.f csyr2k.f ctrmm.f ctrsm.f
88-
chemm.f cherk.f cher2k.f cgemmtr.f)
88+
chemm.f cherk.f cher2k.f cgemmtr.f crotc.f90 scrotc.f90)
8989

90-
set(DBLAS3 dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.f dgemmtr.f)
90+
set(DBLAS3 dgemm.f dsymm.f dsyrk.f dsyr2k.f dtrmm.f dtrsm.f dgemmtr.f drotc.f90)
9191

9292
set(ZBLAS3 zgemm.f zsymm.f zsyrk.f zsyr2k.f ztrmm.f ztrsm.f
93-
zhemm.f zherk.f zher2k.f zgemmtr.f)
93+
zhemm.f zherk.f zher2k.f zgemmtr.f zrotc.f90 dzrotc.f90)
9494

9595

9696
set(SOURCES)

‎BLAS/SRC/crotc.f90‎

Lines changed: 256 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,256 @@
1+
!> \brief \b CROTC applies a chain of rotation sequences to a matrix.
2+
!
3+
! =========== DOCUMENTATION ===========
4+
!
5+
! Online html documentation available at
6+
! http://www.netlib.org/lapack/explore-html/
7+
!
8+
! Definition:
9+
! ===========
10+
!
11+
! subroutine crotc(side, dir, startup, shutdown, m, n, k,&
12+
! A, lda, C, ldc, S, lds)
13+
! .. Scalar Arguments ..
14+
! integer, intent(in) :: m, n, k
15+
! ...
16+
!
17+
!> \par Purpose:
18+
! =============
19+
!>
20+
!> \verbatim
21+
!>
22+
!> CROTC applies a chain of k rotation sequences of length n to a matrix A.
23+
!>
24+
!> Each rotation is specified by a cosine and a sine, stored in the
25+
!> matrices C and S respectively. Rotation G(i,j) is formed by
26+
!> C(i,j) and S(i,j).
27+
!>
28+
!> If side = 'L', rotation G(i,j) is applied to rows i and i+1 of A.
29+
!> [ A(i,j) ] = [ C(i,j) S(i,j) ] [ A(i,j) ]
30+
!> [ A(i+1,j) ] [ -conj(S(i,j)) C(i,j) ] [ A(i+1,j) ]
31+
!> If side = 'R', rotation G(i,j) is applied to columns j and j+1 of A.
32+
!> [ A(i,j) A(i,j+1) ] = [ A(i,j) A(i,j+1) ] [ C(i,j) -conj(S(i,j)) ]
33+
!> [ A(i+1,j) A(i+1,j+1) ] [ A(i+1,j) A(i+1,j+1) ] [ S(i,j) C(i,j) ]
34+
!>
35+
!> \endverbatim
36+
!
37+
! Arguments:
38+
! ==========
39+
!
40+
!> \param[in] side
41+
!> \verbatim
42+
!> side is CHARACTER*1
43+
!> If side = 'L', the rotations are applied to A from the left.
44+
!> If side = 'R', the rotations are applied to A from the right.
45+
!> \endverbatim
46+
!>
47+
!> \param[in] dir
48+
!> \verbatim
49+
!> dir is CHARACTER*1
50+
!> If dir = 'F', the rotations are applied in sequence from the
51+
!> first column/row to the last column/row.
52+
!> If dir = 'B', the rotations are applied in sequence from the
53+
!> last column/row to the first column/row.
54+
!> \endverbatim
55+
!>
56+
!> \param[in] startup
57+
!> \verbatim
58+
!> startup is LOGICAL
59+
!> If startup = .FALSE., the first (k-1) x (k-1) triangle
60+
!> of rotations is not applied.
61+
!> \endverbatim
62+
!>
63+
!> \param[in] shutdown
64+
!> \verbatim
65+
!> shutdown is LOGICAL
66+
!> If shutdown = .FALSE., the last (k-1) x (k-1) triangle
67+
!> of rotations is not applied.
68+
!> \endverbatim
69+
!>
70+
!> \param[in] m
71+
!> \verbatim
72+
!> m is INTEGER
73+
!> If side = 'L', m is the number of columns of A.
74+
!> If side = 'R', m is the number of rows of A.
75+
!> \endverbatim
76+
!>
77+
!> \param[in] n
78+
!> \verbatim
79+
!> n is INTEGER
80+
!> The number of rotations in one sequence.
81+
!> \endverbatim
82+
!>
83+
!> \param[in] k
84+
!> \verbatim
85+
!> k is INTEGER
86+
!> The number of sequences of rotations.
87+
!> \endverbatim
88+
!>
89+
!> \param[in,out] A
90+
!> \verbatim
91+
!> A is COMPLEX array
92+
!> If side = 'L', A has dimension (n+1,m).
93+
!> If side = 'R', A has dimension (m,n+1).
94+
!> The matrix to which the rotations are applied.
95+
!> \endverbatim
96+
!>
97+
!> \param[in] lda
98+
!> \verbatim
99+
!> lda is INTEGER
100+
!> The leading dimension of A.
101+
!> If side = 'L', lda >= n+1.
102+
!> If side = 'R', lda >= m.
103+
!> \endverbatim
104+
!>
105+
!> \param[in,out] C
106+
!> \verbatim
107+
!> C is REAL array, dimension (ldc,k)
108+
!> The matrix containing the cosines of the rotations.
109+
!> \endverbatim
110+
!>
111+
!> \param[in] ldc
112+
!> \verbatim
113+
!> ldc is INTEGER
114+
!> The leading dimension of C.
115+
!> ldc >= n.
116+
!> \endverbatim
117+
!>
118+
!> \param[in,out] S
119+
!> \verbatim
120+
!> S is COMPLEX array, dimension (lds,k)
121+
!> The matrix containing the sines of the rotations.
122+
!> \endverbatim
123+
!>
124+
!> \param[in] lds
125+
!> \verbatim
126+
!> lds is INTEGER
127+
!> The leading dimension of S.
128+
!> lds >= n.
129+
!> \endverbatim
130+
!
131+
! Authors:
132+
! ========
133+
!
134+
!> \author Thijs Steel, KU Leuven, Belgium
135+
!
136+
!> \date October 2024
137+
!
138+
!> \ingroup rotc
139+
!
140+
subroutine crotc(side, dir, startup, shutdown, m, n, k,&
141+
A, lda, C, ldc, S, lds)
142+
! .. Scalar Arguments ..
143+
integer, intent(in) :: m, n, k, lda, ldc, lds
144+
character, intent(in) :: dir, side
145+
logical, intent(in) :: startup, shutdown
146+
! .. Array Arguments ..
147+
complex, intent(inout) :: A(lda,*)
148+
complex, intent(in) :: S(lds,*)
149+
real, intent(in) :: C(ldc,*)
150+
! .. Local Scalars ..
151+
integer i, j, l, j1, j2, incj, incj1, incj2, info
152+
complex temp, sn
153+
real cs
154+
! .. Executable Statements ..
155+
156+
! Test the input parameters
157+
info = 0
158+
if(.not. (side .eq. 'L' .or. side .eq. 'R')) then
159+
info = 1
160+
end if
161+
if(.not. (dir .eq. 'F' .or. dir .eq. 'B')) then
162+
info = 2
163+
end if
164+
if(m .lt. 0) then
165+
info = 5
166+
end if
167+
if(n .lt. 0) then
168+
info = 6
169+
end if
170+
if(k .lt. 0) then
171+
info = 7
172+
end if
173+
if(side .eq. 'L') then
174+
if(lda .lt. n+1) then
175+
info = 9
176+
end if
177+
else
178+
if(lda .lt. m) then
179+
info = 9
180+
end if
181+
end if
182+
if(ldc .lt. n) then
183+
info = 11
184+
end if
185+
if(lds .lt. n) then
186+
info = 13
187+
end if
188+
189+
if(info .ne. 0) then
190+
call xerbla('CROTC ', info)
191+
return
192+
end if
193+
194+
! Determine ranges for loops around C and S
195+
! The range for sequence l is:
196+
! j1+(l-1)*incj1:incj:j2+(l-1)*incj2
197+
if( dir .eq. 'F') then
198+
incj = 1
199+
if(startup) then
200+
j1 = 1
201+
incj1 = 0
202+
else
203+
j1 = k
204+
incj1 = -1
205+
end if
206+
j2 = n
207+
if(shutdown) then
208+
incj2 = 0
209+
else
210+
incj2 = -1
211+
end if
212+
else
213+
incj = -1
214+
j1 = 1
215+
if(startup) then
216+
incj1 = 1
217+
else
218+
incj1 = 0
219+
end if
220+
if(shutdown) then
221+
j2 = 0
222+
incj2 = 0
223+
else
224+
j2 = n-k+1
225+
incj2 = 1
226+
end if
227+
end if
228+
229+
! Apply the rotations
230+
if(side .eq. 'L') then
231+
do l = 1, k
232+
do j = j1+(l-1)*incj1, j2+(l-1)*incj2, incj
233+
cs = C(j,l)
234+
sn = S(j,l)
235+
do i = 1, m
236+
temp = cs*A(i,j) + sn*A(i,j+1)
237+
A(i,j+1) = -conj(sn*A(i,j)) + cs*A(i,j+1)
238+
A(i,j) = temp
239+
end do
240+
end do
241+
end do
242+
else
243+
do l = 1, k
244+
do j = j1+(l-1)*incj1, j2+(l-1)*incj2, incj
245+
cs = C(l,j)
246+
sn = S(l,j)
247+
do i = 1, m
248+
temp = cs*A(j,i) + sn*A(j+1,i)
249+
A(j+1,i) = -conj(sn*A(j,i)) + cs*A(j+1,i)
250+
A(j,i) = temp
251+
end do
252+
end do
253+
end do
254+
end if
255+
256+
end subroutine crotc

0 commit comments

Comments
 (0)