Skip to content

Commit 93515c2

Browse files
martin-frbgmartin-frbg
authored
Merge pull request #5736 from martin-frbg/lapack1221
Follow-up on ?GESVDQ updates from PR1146 (Reference-LAPACK PR 1221)
2 parents d27e98c + 7dde52d commit 93515c2

5 files changed

Lines changed: 405 additions & 199 deletions

File tree

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

Lines changed: 80 additions & 39 deletions
Original file line numberDiff line numberDiff line change
@@ -50,12 +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, 's' ) ) ? m : 1;
55-
lapack_int ncols_u = LAPACKE_lsame( jobu, 'a' ) ? m :
56-
(LAPACKE_lsame( jobu, 's' ) ? MIN(m,n) : 1);
57-
lapack_int nrows_v = LAPACKE_lsame( jobv, 'a' ) ? n :
58-
( LAPACKE_lsame( jobv, 's' ) ? MIN(m,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+
5989
lapack_int lda_t = MAX(1,m);
6090
lapack_int ldu_t = MAX(1,nrows_u);
6191
lapack_int ldv_t = MAX(1,nrows_v);
@@ -73,69 +103,80 @@ lapack_int LAPACKE_cgesvdq_work( int matrix_layout, char joba, char jobp,
73103
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );
74104
return info;
75105
}
76-
if( ldv < n ) {
106+
if( ldv < ncols_v ) {
77107
info = -14;
78108
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );
79109
return info;
80110
}
81111
/* Query optimal working array(s) size if requested */
82-
if( lcwork == -1 ) {
112+
if ( ( liwork == -1 ) || ( lcwork == -1 ) || ( lrwork == -1 ) ) {
83113
LAPACK_cgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
84114
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
85115
cwork, &lcwork, rwork, &lrwork, &info );
86116
return (info < 0) ? (info - 1) : info;
87117
}
88118
/* Allocate memory for temporary array(s) */
89-
a_t = (lapack_complex_float*)LAPACKE_malloc( sizeof(lapack_complex_float) * lda_t * MAX(1,n) );
90-
if( a_t == NULL ) {
91-
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
92-
goto exit_level_0;
93-
}
94-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
95-
u_t = (lapack_complex_float*)
96-
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) );
97127
if( u_t == NULL ) {
98128
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
99129
goto exit_level_1;
100130
}
101131
}
102-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
103-
v_t = (lapack_complex_float*)
104-
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) );
105134
if( v_t == NULL ) {
106135
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
107136
goto exit_level_2;
108137
}
109-
}
138+
110139
/* Transpose input matrices */
111-
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+
112144
/* Call LAPACK function and adjust info */
113-
LAPACK_cgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
114-
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
115-
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 );
116148
if( info < 0 ) {
117149
info = info - 1;
118150
}
151+
119152
/* Transpose output matrices */
120-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
121-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
122-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
123-
u, ldu );
124-
}
125-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
126-
LAPACKE_cge_trans( LAPACK_COL_MAJOR, nrows_v, n, v_t, ldv_t, v,
127-
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+
}
128172
}
173+
129174
/* Release memory and exit */
130-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
131-
LAPACKE_free( v_t );
132-
}
175+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( v_t ); v_t = NULL; }
133176
exit_level_2:
134-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
135-
LAPACKE_free( u_t );
136-
}
177+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( u_t ); u_t = NULL; }
137178
exit_level_1:
138-
LAPACKE_free( a_t );
179+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( a_t ); a_t = NULL; }
139180
exit_level_0:
140181
if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
141182
LAPACKE_xerbla( "LAPACKE_cgesvdq_work", info );

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

Lines changed: 79 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -50,12 +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, 's' ) ) ? m : 1;
55-
lapack_int ncols_u = LAPACKE_lsame( jobu, 'a' ) ? m :
56-
(LAPACKE_lsame( jobu, 's' ) ? MIN(m,n) : 1);
57-
lapack_int nrows_v = LAPACKE_lsame( jobv, 'a' ) ? n :
58-
( LAPACKE_lsame( jobv, 's' ) ? MIN(m,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+
5989
lapack_int lda_t = MAX(1,m);
6090
lapack_int ldu_t = MAX(1,nrows_u);
6191
lapack_int ldv_t = MAX(1,nrows_v);
@@ -73,69 +103,82 @@ lapack_int LAPACKE_dgesvdq_work( int matrix_layout, char joba, char jobp,
73103
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );
74104
return info;
75105
}
76-
if( ldv < n ) {
106+
if( ldv < ncols_v ) {
77107
info = -14;
78108
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );
79109
return info;
80110
}
111+
81112
/* Query optimal working array(s) size if requested */
82-
if( lwork == -1 ) {
113+
if ( ( liwork == -1 ) || ( lwork == -1 ) || ( lrwork == -1 ) ) {
83114
LAPACK_dgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
84115
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
85116
work, &lwork, rwork, &lrwork, &info );
86117
return (info < 0) ? (info - 1) : info;
87118
}
119+
88120
/* Allocate memory for temporary array(s) */
89-
a_t = (double*)LAPACKE_malloc( sizeof(double) * lda_t * MAX(1,n) );
90-
if( a_t == NULL ) {
91-
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
92-
goto exit_level_0;
93-
}
94-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
95-
u_t = (double*)
96-
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) );
97129
if( u_t == NULL ) {
98130
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
99131
goto exit_level_1;
100132
}
101-
}
102-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
103-
v_t = (double*)
104-
LAPACKE_malloc( sizeof(double) * ldv_t * MAX(1,n) );
133+
134+
v_t = (double*)LAPACKE_malloc( sizeof(double) * ldv_t * MAX(1,ncols_v) );
105135
if( v_t == NULL ) {
106136
info = LAPACK_TRANSPOSE_MEMORY_ERROR;
107137
goto exit_level_2;
108138
}
109139
}
140+
110141
/* Transpose input matrices */
111-
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+
112146
/* Call LAPACK function and adjust info */
113147
LAPACK_dgesvdq( &joba, &jobp, &jobr, &jobu, &jobv, &m, &n, a, &lda_t,
114148
s, u, &ldu_t, v, &ldv_t, numrank, iwork, &liwork,
115149
work, &lwork, rwork, &lrwork, &info );
116150
if( info < 0 ) {
117151
info = info - 1;
118152
}
153+
119154
/* Transpose output matrices */
120-
LAPACKE_dge_trans( LAPACK_COL_MAJOR, m, n, a_t, lda_t, a, lda );
121-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
122-
LAPACKE_dge_trans( LAPACK_COL_MAJOR, nrows_u, ncols_u, u_t, ldu_t,
123-
u, ldu );
124-
}
125-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
126-
LAPACKE_dge_trans( LAPACK_COL_MAJOR, nrows_v, n, v_t, ldv_t, v,
127-
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+
}
128174
}
175+
129176
/* Release memory and exit */
130-
if( LAPACKE_lsame( jobv, 'a' ) || LAPACKE_lsame( jobv, 's' ) ) {
131-
LAPACKE_free( v_t );
132-
}
177+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( v_t ); v_t = NULL; }
133178
exit_level_2:
134-
if( LAPACKE_lsame( jobu, 'a' ) || LAPACKE_lsame( jobu, 's' ) ) {
135-
LAPACKE_free( u_t );
136-
}
179+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( u_t ); u_t = NULL; }
137180
exit_level_1:
138-
LAPACKE_free( a_t );
181+
if ( ( m > 0 ) && ( n > 0 ) ) { LAPACKE_free( a_t ); a_t = NULL; }
139182
exit_level_0:
140183
if( info == LAPACK_TRANSPOSE_MEMORY_ERROR ) {
141184
LAPACKE_xerbla( "LAPACKE_dgesvdq_work", info );

0 commit comments

Comments
 (0)