1 // Copyright (C) 2010 Davis E. King (davis@dlib.net) 2 // License: Boost Software License See LICENSE.txt for the full license. 3 #ifndef DLIB_LAPACk_ORMQR_Hh_ 4 #define DLIB_LAPACk_ORMQR_Hh_ 5 6 #include "fortran_id.h" 7 #include "../matrix.h" 8 9 namespace dlib 10 { 11 namespace lapack 12 { 13 namespace binding 14 { 15 extern "C" 16 { 17 void DLIB_FORTRAN_ID(dormqr) (const char *side, const char *trans, const integer *m, const integer *n, 18 const integer *k, const double *a, const integer *lda, const double *tau, 19 double * c_, const integer *ldc, double *work, const integer *lwork, 20 integer *info); 21 22 void DLIB_FORTRAN_ID(sormqr) (const char *side, const char *trans, const integer *m, const integer *n, 23 const integer *k, const float *a, const integer *lda, const float *tau, 24 float * c_, const integer *ldc, float *work, const integer *lwork, 25 integer *info); 26 27 } 28 ormqr(char side,char trans,integer m,integer n,integer k,const double * a,integer lda,const double * tau,double * c_,integer ldc,double * work,integer lwork)29 inline int ormqr (char side, char trans, integer m, integer n, 30 integer k, const double *a, integer lda, const double *tau, 31 double *c_, integer ldc, double *work, integer lwork) 32 { 33 integer info = 0; 34 DLIB_FORTRAN_ID(dormqr)(&side, &trans, &m, &n, 35 &k, a, &lda, tau, 36 c_, &ldc, work, &lwork, &info); 37 return info; 38 } 39 ormqr(char side,char trans,integer m,integer n,integer k,const float * a,integer lda,const float * tau,float * c_,integer ldc,float * work,integer lwork)40 inline int ormqr (char side, char trans, integer m, integer n, 41 integer k, const float *a, integer lda, const float *tau, 42 float *c_, integer ldc, float *work, integer lwork) 43 { 44 integer info = 0; 45 DLIB_FORTRAN_ID(sormqr)(&side, &trans, &m, &n, 46 &k, a, &lda, tau, 47 c_, &ldc, work, &lwork, &info); 48 return info; 49 } 50 51 52 53 } 54 55 // ------------------------------------------------------------------------------------ 56 57 /* -- LAPACK routine (version 3.1) -- */ 58 /* Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */ 59 /* November 2006 */ 60 61 /* .. Scalar Arguments .. */ 62 /* .. */ 63 /* .. Array Arguments .. */ 64 /* .. */ 65 66 /* Purpose */ 67 /* ======= */ 68 69 /* DORMQR overwrites the general real M-by-N matrix C with */ 70 71 /* SIDE = 'L' SIDE = 'R' */ 72 /* TRANS = 'N': Q * C C * Q */ 73 /* TRANS = 'T': Q**T * C C * Q**T */ 74 75 /* where Q is a real orthogonal matrix defined as the product of k */ 76 /* elementary reflectors */ 77 78 /* Q = H(1) H(2) . . . H(k) */ 79 80 /* as returned by DGEQRF. Q is of order M if SIDE = 'L' and of order N */ 81 /* if SIDE = 'R'. */ 82 83 /* Arguments */ 84 /* ========= */ 85 86 /* SIDE (input) CHARACTER*1 */ 87 /* = 'L': apply Q or Q**T from the Left; */ 88 /* = 'R': apply Q or Q**T from the Right. */ 89 90 /* TRANS (input) CHARACTER*1 */ 91 /* = 'N': No transpose, apply Q; */ 92 /* = 'T': Transpose, apply Q**T. */ 93 94 /* M (input) INTEGER */ 95 /* The number of rows of the matrix C. M >= 0. */ 96 97 /* N (input) INTEGER */ 98 /* The number of columns of the matrix C. N >= 0. */ 99 100 /* K (input) INTEGER */ 101 /* The number of elementary reflectors whose product defines */ 102 /* the matrix Q. */ 103 /* If SIDE = 'L', M >= K >= 0; */ 104 /* if SIDE = 'R', N >= K >= 0. */ 105 106 /* A (input) DOUBLE PRECISION array, dimension (LDA,K) */ 107 /* The i-th column must contain the vector which defines the */ 108 /* elementary reflector H(i), for i = 1,2,...,k, as returned by */ 109 /* DGEQRF in the first k columns of its array argument A. */ 110 /* A is modified by the routine but restored on exit. */ 111 112 /* LDA (input) INTEGER */ 113 /* The leading dimension of the array A. */ 114 /* If SIDE = 'L', LDA >= max(1,M); */ 115 /* if SIDE = 'R', LDA >= max(1,N). */ 116 117 /* TAU (input) DOUBLE PRECISION array, dimension (K) */ 118 /* TAU(i) must contain the scalar factor of the elementary */ 119 /* reflector H(i), as returned by DGEQRF. */ 120 121 /* C (input/output) DOUBLE PRECISION array, dimension (LDC,N) */ 122 /* On entry, the M-by-N matrix C. */ 123 /* On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q. */ 124 125 /* LDC (input) INTEGER */ 126 /* The leading dimension of the array C. LDC >= max(1,M). */ 127 128 /* WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) */ 129 /* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. */ 130 131 /* LWORK (input) INTEGER */ 132 /* The dimension of the array WORK. */ 133 /* If SIDE = 'L', LWORK >= max(1,N); */ 134 /* if SIDE = 'R', LWORK >= max(1,M). */ 135 /* For optimum performance LWORK >= N*NB if SIDE = 'L', and */ 136 /* LWORK >= M*NB if SIDE = 'R', where NB is the optimal */ 137 /* blocksize. */ 138 139 /* If LWORK = -1, then a workspace query is assumed; the routine */ 140 /* only calculates the optimal size of the WORK array, returns */ 141 /* this value as the first entry of the WORK array, and no error */ 142 /* message related to LWORK is issued by XERBLA. */ 143 144 /* INFO (output) INTEGER */ 145 /* = 0: successful exit */ 146 /* < 0: if INFO = -i, the i-th argument had an illegal value */ 147 148 // ------------------------------------------------------------------------------------ 149 150 template < 151 typename T, 152 long NR1, long NR2, long NR3, 153 long NC1, long NC2, long NC3, 154 typename MM, 155 typename C_LAYOUT 156 > ormqr(char side,char trans,const matrix<T,NR1,NC1,MM,column_major_layout> & a,const matrix<T,NR2,NC2,MM,column_major_layout> & tau,matrix<T,NR3,NC3,MM,C_LAYOUT> & c)157 int ormqr ( 158 char side, 159 char trans, 160 const matrix<T,NR1,NC1,MM,column_major_layout>& a, 161 const matrix<T,NR2,NC2,MM,column_major_layout>& tau, 162 matrix<T,NR3,NC3,MM,C_LAYOUT>& c 163 ) 164 { 165 long m = c.nr(); 166 long n = c.nc(); 167 const long k = a.nc(); 168 long ldc; 169 if (is_same_type<C_LAYOUT,column_major_layout>::value) 170 { 171 ldc = c.nr(); 172 } 173 else 174 { 175 // Since lapack expects c to be in column major layout we have to 176 // do something to make this work. Since a row major layout matrix 177 // will look just like a transposed C we can just swap a few things around. 178 179 ldc = c.nc(); 180 swap(m,n); 181 182 if (side == 'L') 183 side = 'R'; 184 else 185 side = 'L'; 186 187 if (trans == 'T') 188 trans = 'N'; 189 else 190 trans = 'T'; 191 } 192 193 matrix<T,0,1,MM,column_major_layout> work; 194 195 // figure out how big the workspace needs to be. 196 T work_size = 1; 197 int info = binding::ormqr(side, trans, m, n, 198 k, &a(0,0), a.nr(), &tau(0,0), 199 &c(0,0), ldc, &work_size, -1); 200 201 if (info != 0) 202 return info; 203 204 if (work.size() < work_size) 205 work.set_size(static_cast<long>(work_size), 1); 206 207 // compute the actual result 208 info = binding::ormqr(side, trans, m, n, 209 k, &a(0,0), a.nr(), &tau(0,0), 210 &c(0,0), ldc, &work(0,0), work.size()); 211 212 return info; 213 } 214 215 // ------------------------------------------------------------------------------------ 216 217 } 218 219 } 220 221 // ---------------------------------------------------------------------------------------- 222 223 #endif // DLIB_LAPACk_ORMQR_Hh_ 224 225