Skip to content

Commit 7dde52d

Browse files
martin-frbgmartin-frbg
authored
Follow-up on ?GESVDQ updates from PR 1146 (Reference-LAPACK PR 1221)
1 parent c6e4d17 commit 7dde52d

4 files changed

Lines changed: 318 additions & 218 deletions

File tree

‎lapack-netlib/LAPACKE/src/lapacke_cgesvdq_work.c‎

Lines changed: 80 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -50,17 +50,42 @@ lapack_int LAPACKE_cgesvdq_work( int matrix_layout, char joba, char jobp,
5050
info = info - 1;
5151
}
5252
} else if( matrix_layout == LAPACK_ROW_MAJOR ) {
53-
lapack_int nrows_u = LAPACKE_lsame( jobu, 'a' ) ||
54-
LAPACKE_lsame( jobu, 'u' ) ||
55-
LAPACKE_lsame( jobu, 'r' ) ||
56-
LAPACKE_lsame( jobu, 'f' ) ||
57-
LAPACKE_lsame( jobu, 's' ) ? m : 1;
58-
lapack_int ncols_u = LAPACKE_lsame( jobu, 'a' ) ? m :
59-
( LAPACKE_lsame( jobu, 's' ) ||
60-
(LAPACKE_lsame( jobu, 'u' ) ) ? MIN(m,n) : 1);
61-
lapack_int nrows_v = ( LAPACKE_lsame( jobv, 'a' ) ||
62-
LAPACKE_lsame( jobv, 'v' ) ||
63-
LAPACKE_lsame( jobv, 'r' )) ? n : 1;
53+
54+
lapack_int nrows_u;
55+
lapack_int ncols_u;
56+
lapack_int nrows_v;
57+
lapack_int ncols_v;
58+
59+
if( LAPACKE_lsame( jobu, 'a' ) ) {
60+
nrows_u = m;
61+
ncols_u = m;
62+
}
63+
else if( LAPACKE_lsame( jobu, 's' ) ||
64+
LAPACKE_lsame( jobu, 'u' ) ||
65+
LAPACKE_lsame( jobu, 'r' ) ) {
66+
nrows_u = m;
67+
ncols_u = n;
68+
}
69+
else if( LAPACKE_lsame( jobu, 'f' ) ) {
70+
nrows_u = n;
71+
ncols_u = n;
72+
} else {
73+
nrows_u = 1;
74+
ncols_u = 1;
75+
}
76+
77+
/* in the case joba == 'e', v_t is used as a workspace */
78+
if( LAPACKE_lsame( jobv, 'a' ) ||
79+
LAPACKE_lsame( jobv, 'v' ) ||
80+
LAPACKE_lsame( jobv, 'r' ) ||
81+
LAPACKE_lsame( joba, 'e' ) ) {
82+
nrows_v = n;
83+
ncols_v = n;
84+
} else {
85+
nrows_v = 1;
86+
ncols_v = 1;
87+
}
88+
6489
lapack_int lda_t = MAX(1,m);
6590
lapack_int ldu_t = MAX(1,nrows_u);
6691
lapack_int ldv_t = MAX(1,nrows_v);
@@ -78,81 +103,80 @@ lapack_int LAPACKE_cgesvdq_work( int matrix_layout, char joba, char jobp,
78103
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );
79104
return info;
80105
}
81-
if( ldv < n ) {
106+
if( ldv < ncols_v ) {
82107
info = -14;
83108
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );
84109
return info;
85110
}
86111
/* Query optimal working array(s) size if requested */
87-
if( lcwork == -1 ) {
112+
if ( ( liwork == -1 ) || ( lcwork == -1 ) || ( lrwork == -1 ) ) {
88113
LAPACK_cgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
89114
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
90115
cwork, &lcwork, rwork, &lrwork, &info );
91116
return (info < 0) ? (info - 1) : info;
92117
}
93118
/* Allocate memory for temporary array(s) */
94-
a_t = (lapack_complex_float*)LAPACKE_malloc( sizeof(lapack_complex_float) * lda_t * MAX(1,n) );
95-
if( a_t == NULL ) {
96-
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
97-
goto exit_level_0;
98-
}
99-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ||
100-
LAPACKE_lsame( jobu, 'u' ) ||
101-
LAPACKE_lsame( jobu, 'r' ) ||
102-
LAPACKE_lsame( jobu, 'f' ) ) {
103-
u_t = (lapack_complex_float*)
104-
LAPACKE_malloc( sizeof(lapack_complex_float) * ldu_t * MAX(1,ncols_u) );
119+
if ( ( m > 0 ) && ( n > 0 ) ){
120+
a_t = (lapack_complex_float*)LAPACKE_malloc( sizeof(lapack_complex_float) * lda_t * MAX(1,n) );
121+
if( a_t == NULL ) {
122+
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
123+
goto exit_level_0;
124+
}
125+
126+
u_t = (lapack_complex_float*)LAPACKE_malloc( sizeof(lapack_complex_float) * ldu_t * MAX(1,ncols_u) );
105127
if( u_t == NULL ) {
106128
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
107129
goto exit_level_1;
108130
}
109131
}
110-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
111-
LAPACKE_lsame( jobv, 'r' ) ) {
112-
v_t = (lapack_complex_float*)
113-
LAPACKE_malloc( sizeof(lapack_complex_float) * ldv_t * MAX(1,n) );
132+
133+
v_t = (lapack_complex_float*)LAPACKE_malloc( sizeof(lapack_complex_float) * ldv_t * MAX(1,ncols_v) );
114134
if( v_t == NULL ) {
115135
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
116136
goto exit_level_2;
117137
}
118-
}
138+
119139
/* Transpose input matrices */
120-
LAPACKE_cge_trans( matrix_layout, m, n, a, lda, a_t, lda_t );
140+
if ( ( m > 0 ) && ( n > 0 ) ){
141+
LAPACKE_cge_trans( matrix_layout, m, n, a, lda, a_t, lda_t );
142+
}
143+
121144
/* Call LAPACK function and adjust info */
122-
LAPACK_cgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
123-
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
124-
cwork, &lcwork, rwork, &lrwork, &info );
145+
LAPACK_cgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
146+
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
147+
cwork, &lcwork, rwork, &lrwork, &info );
125148
if( info < 0 ) {
126149
info = info - 1;
127150
}
151+
128152
/* Transpose output matrices */
129-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
130-
if( LAPACKE_lsame( jobu, 'a' ) ||LAPACKE_lsame( jobu, 's' ) ||
131-
LAPACKE_lsame( jobu, 'u' ) ||
132-
LAPACKE_lsame( jobu, 'r' ) ||
133-
LAPACKE_lsame( jobu, 'f' ) ) {
134-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
135-
u, ldu );
136-
}
137-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
138-
LAPACKE_lsame( jobv, 'r' )) {
139-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_v, n, v_t, ldv_t, v,
140-
ldv );
153+
if ( ( m > 0 ) && ( n > 0 ) ){
154+
LAPACKE_cge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
155+
156+
if( LAPACKE_lsame( jobu, 'a' ) ||
157+
LAPACKE_lsame( jobu, 's' ) ||
158+
LAPACKE_lsame( jobu, 'u' ) ||
159+
LAPACKE_lsame( jobu, 'r' ) ||
160+
LAPACKE_lsame( jobu, 'f' ) ) {
161+
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
162+
u, ldu );
163+
}
164+
165+
/* we do not transpose v_t back to v in the case (joba == 'e') because, in this case, v_t is used as a workspace */
166+
if( LAPACKE_lsame( jobv, 'a' ) ||
167+
LAPACKE_lsame( jobv, 'v' ) ||
168+
LAPACKE_lsame( jobv, 'r' )) {
169+
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_v, ncols_v, v_t, ldv_t, v,
170+
ldv );
171+
}
141172
}
173+
142174
/* Release memory and exit */
143-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
144-
LAPACKE_lsame( jobv, 'r' ) ) {
145-
LAPACKE_free( v_t );
146-
}
175+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( v_t ); v_t = NULL; }
147176
exit_level_2:
148-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ||
149-
LAPACKE_lsame( jobu, 'u' ) ||
150-
LAPACKE_lsame( jobu, 'r' ) ||
151-
LAPACKE_lsame( jobu, 'f' ) ) {
152-
LAPACKE_free( u_t );
153-
}
177+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( u_t ); u_t = NULL; }
154178
exit_level_1:
155-
LAPACKE_free( a_t );
179+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( a_t ); a_t = NULL; }
156180
exit_level_0:
157181
if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
158182
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );

‎lapack-netlib/LAPACKE/src/lapacke_dgesvdq_work.c‎

Lines changed: 79 additions & 53 deletions
Original file line numberDiff line numberDiff line change
@@ -50,17 +50,42 @@ lapack_int LAPACKE_dgesvdq_work( int matrix_layout, char joba, char jobp,
5050
info = info - 1;
5151
}
5252
} else if( matrix_layout == LAPACK_ROW_MAJOR ) {
53-
lapack_int nrows_u = ( LAPACKE_lsame( jobu, 'a' ) ||
54-
LAPACKE_lsame( jobu, 'u' ) ||
55-
LAPACKE_lsame( jobu, 'r' ) ||
56-
LAPACKE_lsame( jobu, 'f' ) ||
57-
LAPACKE_lsame( jobu, 's' ) ) ? m : 1;
58-
lapack_int ncols_u = LAPACKE_lsame( jobu, 'a' ) ? m :
59-
( (LAPACKE_lsame( jobu, 's' ) ||
60-
LAPACKE_lsame( jobu, 'u' ) ) ? MIN(m,n) : 1);
61-
lapack_int nrows_v = ( LAPACKE_lsame( jobv, 'a' ) ||
62-
LAPACKE_lsame( jobu, 'v' ) ||
63-
LAPACKE_lsame( jobu, 'r' )) ? n : 1;
53+
54+
lapack_int nrows_u;
55+
lapack_int ncols_u;
56+
lapack_int nrows_v;
57+
lapack_int ncols_v;
58+
59+
if( LAPACKE_lsame( jobu, 'a' ) ) {
60+
nrows_u = m;
61+
ncols_u = m;
62+
}
63+
else if( LAPACKE_lsame( jobu, 's' ) ||
64+
LAPACKE_lsame( jobu, 'u' ) ||
65+
LAPACKE_lsame( jobu, 'r' ) ) {
66+
nrows_u = m;
67+
ncols_u = n;
68+
}
69+
else if( LAPACKE_lsame( jobu, 'f' ) ) {
70+
nrows_u = n;
71+
ncols_u = n;
72+
} else {
73+
nrows_u = 1;
74+
ncols_u = 1;
75+
}
76+
77+
/* in the case joba == 'e', v_t is used as a workspace */
78+
if( LAPACKE_lsame( jobv, 'a' ) ||
79+
LAPACKE_lsame( jobv, 'v' ) ||
80+
LAPACKE_lsame( jobv, 'r' ) ||
81+
LAPACKE_lsame( joba, 'e' ) ) {
82+
nrows_v = n;
83+
ncols_v = n;
84+
} else {
85+
nrows_v = 1;
86+
ncols_v = 1;
87+
}
88+
6489
lapack_int lda_t = MAX(1,m);
6590
lapack_int ldu_t = MAX(1,nrows_u);
6691
lapack_int ldv_t = MAX(1,nrows_v);
@@ -78,81 +103,82 @@ lapack_int LAPACKE_dgesvdq_work( int matrix_layout, char joba, char jobp,
78103
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );
79104
return info;
80105
}
81-
if( ldv < n ) {
106+
if( ldv < ncols_v ) {
82107
info = -14;
83108
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );
84109
return info;
85110
}
111+
86112
/* Query optimal working array(s) size if requested */
87-
if( lwork == -1 ) {
113+
if ( ( liwork == -1 ) || ( lwork == -1 ) || ( lrwork == -1 ) ) {
88114
LAPACK_dgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
89115
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
90116
work, &lwork, rwork, &lrwork, &info );
91117
return (info < 0) ? (info - 1) : info;
92118
}
119+
93120
/* Allocate memory for temporary array(s) */
94-
a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * MAX(1,n) );
95-
if( a_t == NULL ) {
96-
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
97-
goto exit_level_0;
98-
}
99-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ||
100-
LAPACKE_lsame( jobu, 'u' ) ||
101-
LAPACKE_lsame( jobu, 'r' ) ||
102-
LAPACKE_lsame( jobu, 'f' ) ) {
103-
u_t = (double*)
104-
LAPACKE_malloc( sizeof(double) * ldu_t * MAX(1,ncols_u) );
121+
if ( ( m > 0 ) && ( n > 0 ) ){
122+
a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * MAX(1,n) );
123+
if( a_t == NULL ) {
124+
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
125+
goto exit_level_0;
126+
}
127+
128+
u_t = (double*)LAPACKE_malloc( sizeof(double) * ldu_t * MAX(1,ncols_u) );
105129
if( u_t == NULL ) {
106130
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
107131
goto exit_level_1;
108132
}
109-
}
110-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
111-
LAPACKE_lsame( jobu, 'r' ) ) {
112-
v_t = (double*)
113-
LAPACKE_malloc( sizeof(double) * ldv_t * MAX(1,n) );
133+
134+
v_t = (double*)LAPACKE_malloc( sizeof(double) * ldv_t * MAX(1,ncols_v) );
114135
if( v_t == NULL ) {
115136
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
116137
goto exit_level_2;
117138
}
118139
}
140+
119141
/* Transpose input matrices */
120-
LAPACKE_dge_trans( matrix_layout, m, n, a, lda, a_t, lda_t );
142+
if ( ( m > 0 ) && ( n > 0 ) ){
143+
LAPACKE_dge_trans( matrix_layout, m, n, a, lda, a_t, lda_t );
144+
}
145+
121146
/* Call LAPACK function and adjust info */
122147
LAPACK_dgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
123148
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
124149
work, &lwork, rwork, &lrwork, &info );
125150
if( info < 0 ) {
126151
info = info - 1;
127152
}
153+
128154
/* Transpose output matrices */
129-
LAPACKE_dge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
130-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ||
131-
LAPACKE_lsame( jobu, 'u' ) ||
132-
LAPACKE_lsame( jobu, 'r' ) ||
133-
LAPACKE_lsame( jobu, 'f' ) ) {
134-
LAPACKE_dge_trans)( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
135-
u, ldu );
136-
}
137-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
138-
LAPACKE_lsame( jobu, 'r' )) {
139-
LAPACKE_dge_trans( LAPACK_COL_MAJOR, nrows_v, n, v_t, ldv_t, v,
140-
ldv );
155+
if ( ( m > 0 ) && ( n > 0 ) ){
156+
LAPACKE_dge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
157+
158+
if( LAPACKE_lsame( jobu, 'a' ) ||
159+
LAPACKE_lsame( jobu, 's' ) ||
160+
LAPACKE_lsame( jobu, 'u' ) ||
161+
LAPACKE_lsame( jobu, 'r' ) ||
162+
LAPACKE_lsame( jobu, 'f' ) ) {
163+
LAPACKE_dge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
164+
u, ldu );
165+
}
166+
167+
/* we do not transpose v_t back to v in the case (joba == 'e') because, in this case, v_t is used as a workspace */
168+
if( LAPACKE_lsame( jobv, 'a' ) ||
169+
LAPACKE_lsame( jobv, 'v' ) ||
170+
LAPACKE_lsame( jobv, 'r' )) {
171+
LAPACKE_dge_trans( LAPACK_COL_MAJOR, nrows_v, ncols_v, v_t, ldv_t, v,
172+
ldv );
173+
}
141174
}
175+
142176
/* Release memory and exit */
143-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 'v' ) ||
144-
LAPACKE_lsame( jobu, 'r' ) ) {
145-
LAPACKE_free( v_t );
146-
}
177+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( v_t ); v_t = NULL; }
147178
exit_level_2:
148-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ||
149-
LAPACKE_lsame( jobu, 'u' ) ||
150-
LAPACKE_lsame( jobu, 'r' ) ||
151-
LAPACKE_lsame( jobu, 'f' ) ) {
152-
LAPACKE_free( u_t );
153-
}
179+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( u_t ); u_t = NULL; }
154180
exit_level_1:
155-
LAPACKE_free( a_t );
181+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( a_t ); a_t = NULL; }
156182
exit_level_0:
157183
if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
158184
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );

0 commit comments

Comments
 (0)