1 /* -- translated by f2c (version 20191129). 2 You must link the resulting object file with libf2c: 3 on Microsoft Windows system, link with libf2c.lib; 4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm 5 or, if you install libf2c.a in a standard place, with -lf2c -lm 6 -- in that order, at the end of the command line, as in 7 cc *.o -lf2c -lm 8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g., 9 10 http://www.netlib.org/f2c/libf2c.zip 11 */ 12 13 #include "f2c.h" 14 15 /* Table of constant values */ 16 17 static integer c__1 = 1; 18 static integer c_n1 = -1; 19 static integer c__3 = 3; 20 static integer c__2 = 2; 21 static integer c__65 = 65; 22 static doublereal c_b25 = -1.; 23 static doublereal c_b26 = 1.; 24 25 /* > \brief \b DGEHRD 26 27 =========== DOCUMENTATION =========== 28 29 Online html documentation available at 30 http://www.netlib.org/lapack/explore-html/ 31 32 > \htmlonly 33 > Download DGEHRD + dependencies 34 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgehrd. 35 f"> 36 > [TGZ]</a> 37 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgehrd. 38 f"> 39 > [ZIP]</a> 40 > <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgehrd. 41 f"> 42 > [TXT]</a> 43 > \endhtmlonly 44 45 Definition: 46 =========== 47 48 SUBROUTINE DGEHRD( N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO ) 49 50 INTEGER IHI, ILO, INFO, LDA, LWORK, N 51 DOUBLE PRECISION A( LDA, * ), TAU( * ), WORK( * ) 52 53 54 > \par Purpose: 55 ============= 56 > 57 > \verbatim 58 > 59 > DGEHRD reduces a real general matrix A to upper Hessenberg form H by 60 > an orthogonal similarity transformation: Q**T * A * Q = H . 61 > \endverbatim 62 63 Arguments: 64 ========== 65 66 > \param[in] N 67 > \verbatim 68 > N is INTEGER 69 > The order of the matrix A. N >= 0. 70 > \endverbatim 71 > 72 > \param[in] ILO 73 > \verbatim 74 > ILO is INTEGER 75 > \endverbatim 76 > 77 > \param[in] IHI 78 > \verbatim 79 > IHI is INTEGER 80 > 81 > It is assumed that A is already upper triangular in rows 82 > and columns 1:ILO-1 and IHI+1:N. ILO and IHI are normally 83 > set by a previous call to DGEBAL; otherwise they should be 84 > set to 1 and N respectively. See Further Details. 85 > 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 86 > \endverbatim 87 > 88 > \param[in,out] A 89 > \verbatim 90 > A is DOUBLE PRECISION array, dimension (LDA,N) 91 > On entry, the N-by-N general matrix to be reduced. 92 > On exit, the upper triangle and the first subdiagonal of A 93 > are overwritten with the upper Hessenberg matrix H, and the 94 > elements below the first subdiagonal, with the array TAU, 95 > represent the orthogonal matrix Q as a product of elementary 96 > reflectors. See Further Details. 97 > \endverbatim 98 > 99 > \param[in] LDA 100 > \verbatim 101 > LDA is INTEGER 102 > The leading dimension of the array A. LDA >= max(1,N). 103 > \endverbatim 104 > 105 > \param[out] TAU 106 > \verbatim 107 > TAU is DOUBLE PRECISION array, dimension (N-1) 108 > The scalar factors of the elementary reflectors (see Further 109 > Details). Elements 1:ILO-1 and IHI:N-1 of TAU are set to 110 > zero. 111 > \endverbatim 112 > 113 > \param[out] WORK 114 > \verbatim 115 > WORK is DOUBLE PRECISION array, dimension (LWORK) 116 > On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 117 > \endverbatim 118 > 119 > \param[in] LWORK 120 > \verbatim 121 > LWORK is INTEGER 122 > The length of the array WORK. LWORK >= max(1,N). 123 > For optimum performance LWORK >= N*NB, where NB is the 124 > optimal blocksize. 125 > 126 > If LWORK = -1, then a workspace query is assumed; the routine 127 > only calculates the optimal size of the WORK array, returns 128 > this value as the first entry of the WORK array, and no error 129 > message related to LWORK is issued by XERBLA. 130 > \endverbatim 131 > 132 > \param[out] INFO 133 > \verbatim 134 > INFO is INTEGER 135 > = 0: successful exit 136 > < 0: if INFO = -i, the i-th argument had an illegal value. 137 > \endverbatim 138 139 Authors: 140 ======== 141 142 > \author Univ. of Tennessee 143 > \author Univ. of California Berkeley 144 > \author Univ. of Colorado Denver 145 > \author NAG Ltd. 146 147 > \date November 2011 148 149 > \ingroup doubleGEcomputational 150 151 > \par Further Details: 152 ===================== 153 > 154 > \verbatim 155 > 156 > The matrix Q is represented as a product of (ihi-ilo) elementary 157 > reflectors 158 > 159 > Q = H(ilo) H(ilo+1) . . . H(ihi-1). 160 > 161 > Each H(i) has the form 162 > 163 > H(i) = I - tau * v * v**T 164 > 165 > where tau is a real scalar, and v is a real vector with 166 > v(1:i) = 0, v(i+1) = 1 and v(ihi+1:n) = 0; v(i+2:ihi) is stored on 167 > exit in A(i+2:ihi,i), and tau in TAU(i). 168 > 169 > The contents of A are illustrated by the following example, with 170 > n = 7, ilo = 2 and ihi = 6: 171 > 172 > on entry, on exit, 173 > 174 > ( a a a a a a a ) ( a a h h h h a ) 175 > ( a a a a a a ) ( a h h h h a ) 176 > ( a a a a a a ) ( h h h h h h ) 177 > ( a a a a a a ) ( v2 h h h h h ) 178 > ( a a a a a a ) ( v2 v3 h h h h ) 179 > ( a a a a a a ) ( v2 v3 v4 h h h ) 180 > ( a ) ( a ) 181 > 182 > where a denotes an element of the original matrix A, h denotes a 183 > modified element of the upper Hessenberg matrix H, and vi denotes an 184 > element of the vector defining H(i). 185 > 186 > This file is a slight modification of LAPACK-3.0's DGEHRD 187 > subroutine incorporating improvements proposed by Quintana-Orti and 188 > Van de Geijn (2006). (See DLAHR2.) 189 > \endverbatim 190 > 191 ===================================================================== igraphdgehrd_(integer * n,integer * ilo,integer * ihi,doublereal * a,integer * lda,doublereal * tau,doublereal * work,integer * lwork,integer * info)192 Subroutine */ int igraphdgehrd_(integer *n, integer *ilo, integer *ihi, 193 doublereal *a, integer *lda, doublereal *tau, doublereal *work, 194 integer *lwork, integer *info) 195 { 196 /* System generated locals */ 197 integer a_dim1, a_offset, i__1, i__2, i__3, i__4; 198 199 /* Local variables */ 200 integer i__, j; 201 doublereal t[4160] /* was [65][64] */; 202 integer ib; 203 doublereal ei; 204 integer nb, nh, nx, iws; 205 extern /* Subroutine */ int igraphdgemm_(char *, char *, integer *, integer *, 206 integer *, doublereal *, doublereal *, integer *, doublereal *, 207 integer *, doublereal *, doublereal *, integer *); 208 integer nbmin, iinfo; 209 extern /* Subroutine */ int igraphdtrmm_(char *, char *, char *, char *, 210 integer *, integer *, doublereal *, doublereal *, integer *, 211 doublereal *, integer *), igraphdaxpy_( 212 integer *, doublereal *, doublereal *, integer *, doublereal *, 213 integer *), igraphdgehd2_(integer *, integer *, integer *, doublereal *, 214 integer *, doublereal *, doublereal *, integer *), igraphdlahr2_( 215 integer *, integer *, integer *, doublereal *, integer *, 216 doublereal *, doublereal *, integer *, doublereal *, integer *), 217 igraphdlarfb_(char *, char *, char *, char *, integer *, integer *, 218 integer *, doublereal *, integer *, doublereal *, integer *, 219 doublereal *, integer *, doublereal *, integer *), igraphxerbla_(char *, integer *, ftnlen); 220 extern integer igraphilaenv_(integer *, char *, char *, integer *, integer *, 221 integer *, integer *, ftnlen, ftnlen); 222 integer ldwork, lwkopt; 223 logical lquery; 224 225 226 /* -- LAPACK computational routine (version 3.4.0) -- 227 -- LAPACK is a software package provided by Univ. of Tennessee, -- 228 -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 229 November 2011 230 231 232 ===================================================================== 233 234 235 Test the input parameters 236 237 Parameter adjustments */ 238 a_dim1 = *lda; 239 a_offset = 1 + a_dim1; 240 a -= a_offset; 241 --tau; 242 --work; 243 244 /* Function Body */ 245 *info = 0; 246 /* Computing MIN */ 247 i__1 = 64, i__2 = igraphilaenv_(&c__1, "DGEHRD", " ", n, ilo, ihi, &c_n1, ( 248 ftnlen)6, (ftnlen)1); 249 nb = min(i__1,i__2); 250 lwkopt = *n * nb; 251 work[1] = (doublereal) lwkopt; 252 lquery = *lwork == -1; 253 if (*n < 0) { 254 *info = -1; 255 } else if (*ilo < 1 || *ilo > max(1,*n)) { 256 *info = -2; 257 } else if (*ihi < min(*ilo,*n) || *ihi > *n) { 258 *info = -3; 259 } else if (*lda < max(1,*n)) { 260 *info = -5; 261 } else if (*lwork < max(1,*n) && ! lquery) { 262 *info = -8; 263 } 264 if (*info != 0) { 265 i__1 = -(*info); 266 igraphxerbla_("DGEHRD", &i__1, (ftnlen)6); 267 return 0; 268 } else if (lquery) { 269 return 0; 270 } 271 272 /* Set elements 1:ILO-1 and IHI:N-1 of TAU to zero */ 273 274 i__1 = *ilo - 1; 275 for (i__ = 1; i__ <= i__1; ++i__) { 276 tau[i__] = 0.; 277 /* L10: */ 278 } 279 i__1 = *n - 1; 280 for (i__ = max(1,*ihi); i__ <= i__1; ++i__) { 281 tau[i__] = 0.; 282 /* L20: */ 283 } 284 285 /* Quick return if possible */ 286 287 nh = *ihi - *ilo + 1; 288 if (nh <= 1) { 289 work[1] = 1.; 290 return 0; 291 } 292 293 /* Determine the block size 294 295 Computing MIN */ 296 i__1 = 64, i__2 = igraphilaenv_(&c__1, "DGEHRD", " ", n, ilo, ihi, &c_n1, ( 297 ftnlen)6, (ftnlen)1); 298 nb = min(i__1,i__2); 299 nbmin = 2; 300 iws = 1; 301 if (nb > 1 && nb < nh) { 302 303 /* Determine when to cross over from blocked to unblocked code 304 (last block is always handled by unblocked code) 305 306 Computing MAX */ 307 i__1 = nb, i__2 = igraphilaenv_(&c__3, "DGEHRD", " ", n, ilo, ihi, &c_n1, ( 308 ftnlen)6, (ftnlen)1); 309 nx = max(i__1,i__2); 310 if (nx < nh) { 311 312 /* Determine if workspace is large enough for blocked code */ 313 314 iws = *n * nb; 315 if (*lwork < iws) { 316 317 /* Not enough workspace to use optimal NB: determine the 318 minimum value of NB, and reduce NB or force use of 319 unblocked code 320 321 Computing MAX */ 322 i__1 = 2, i__2 = igraphilaenv_(&c__2, "DGEHRD", " ", n, ilo, ihi, & 323 c_n1, (ftnlen)6, (ftnlen)1); 324 nbmin = max(i__1,i__2); 325 if (*lwork >= *n * nbmin) { 326 nb = *lwork / *n; 327 } else { 328 nb = 1; 329 } 330 } 331 } 332 } 333 ldwork = *n; 334 335 if (nb < nbmin || nb >= nh) { 336 337 /* Use unblocked code below */ 338 339 i__ = *ilo; 340 341 } else { 342 343 /* Use blocked code */ 344 345 i__1 = *ihi - 1 - nx; 346 i__2 = nb; 347 for (i__ = *ilo; i__2 < 0 ? i__ >= i__1 : i__ <= i__1; i__ += i__2) { 348 /* Computing MIN */ 349 i__3 = nb, i__4 = *ihi - i__; 350 ib = min(i__3,i__4); 351 352 /* Reduce columns i:i+ib-1 to Hessenberg form, returning the 353 matrices V and T of the block reflector H = I - V*T*V**T 354 which performs the reduction, and also the matrix Y = A*V*T */ 355 356 igraphdlahr2_(ihi, &i__, &ib, &a[i__ * a_dim1 + 1], lda, &tau[i__], t, & 357 c__65, &work[1], &ldwork); 358 359 /* Apply the block reflector H to A(1:ihi,i+ib:ihi) from the 360 right, computing A := A - Y * V**T. V(i+ib,ib-1) must be set 361 to 1 */ 362 363 ei = a[i__ + ib + (i__ + ib - 1) * a_dim1]; 364 a[i__ + ib + (i__ + ib - 1) * a_dim1] = 1.; 365 i__3 = *ihi - i__ - ib + 1; 366 igraphdgemm_("No transpose", "Transpose", ihi, &i__3, &ib, &c_b25, & 367 work[1], &ldwork, &a[i__ + ib + i__ * a_dim1], lda, & 368 c_b26, &a[(i__ + ib) * a_dim1 + 1], lda); 369 a[i__ + ib + (i__ + ib - 1) * a_dim1] = ei; 370 371 /* Apply the block reflector H to A(1:i,i+1:i+ib-1) from the 372 right */ 373 374 i__3 = ib - 1; 375 igraphdtrmm_("Right", "Lower", "Transpose", "Unit", &i__, &i__3, &c_b26, 376 &a[i__ + 1 + i__ * a_dim1], lda, &work[1], &ldwork); 377 i__3 = ib - 2; 378 for (j = 0; j <= i__3; ++j) { 379 igraphdaxpy_(&i__, &c_b25, &work[ldwork * j + 1], &c__1, &a[(i__ + 380 j + 1) * a_dim1 + 1], &c__1); 381 /* L30: */ 382 } 383 384 /* Apply the block reflector H to A(i+1:ihi,i+ib:n) from the 385 left */ 386 387 i__3 = *ihi - i__; 388 i__4 = *n - i__ - ib + 1; 389 igraphdlarfb_("Left", "Transpose", "Forward", "Columnwise", &i__3, & 390 i__4, &ib, &a[i__ + 1 + i__ * a_dim1], lda, t, &c__65, &a[ 391 i__ + 1 + (i__ + ib) * a_dim1], lda, &work[1], &ldwork); 392 /* L40: */ 393 } 394 } 395 396 /* Use unblocked code to reduce the rest of the matrix */ 397 398 igraphdgehd2_(n, &i__, ihi, &a[a_offset], lda, &tau[1], &work[1], &iinfo); 399 work[1] = (doublereal) iws; 400 401 return 0; 402 403 /* End of DGEHRD */ 404 405 } /* igraphdgehrd_ */ 406 407