通过反射因子生成正交矩阵Q,其中反射因子由SYTRD计算得来。
C Interface:
sorgtr_(const char *uplo, const int *n, float *a, const int *lda, const float *tau, float *work, const int *lwork, int *info);
dorgtr_(const char *uplo, const int *n, double *a, const int *lda, const double *tau, double *work, const int *lwork, int *info);
cungtr_(const char *uplo, const int *n, float _Complex *a, const int *lda, const float _Complex *tau, float _Complex *work, const int *lwork, int *info);
zungtr_(const char *uplo, const int *n, double _Complex *a, const int *lda, const double _Complex *tau, double _Complex *work, const int *lwork, int *info);
Fortran Interface:
SORGTR(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO);
DORGTR(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO);
CUNGTR(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO);
ZUNGTR(UPLO, N, A, LDA, TAU, WORK, LWORK, INFO);
参数 |
类型 |
说明 |
输入/输出 |
---|---|---|---|
uplo |
字符型 |
|
输入 |
n |
整数型 |
矩阵A的维度, n≥0。 |
输入 |
a |
|
大小为lda*n。
|
输入,输出 |
lda |
整数型 |
矩阵A的主维,lda≥max(1, n)。 |
输入 |
tau |
|
大小为n-1。 该数组表示H的标量因子。 |
输入 |
work |
|
工作数组,大小为max(1, lwork)。 info=0时,work(0)返回最优lwork大小。 |
输出 |
lwork |
整数型 |
work的大小。
|
输入 |
Info |
整数型 |
|
输出 |
#include "klapack.h"
C Interface:
char uplo = 'L'; int n = 4; int lda = n; int info = 0; double tau[4] = {1.003949, 1.125229, 1.978923, 0.000000}; double *work = NULL; double qwork; int lwork = -1; double *d = (double*)malloc(n * sizeof(double)); double *e = (double*)malloc((n - 1) * sizeof(double)); /* * tau: * 1.003949 1.125229 1.978923 0.000000 * A (4x4, stored in column-major): * 7.027 8.710 1.015 6.929 * 8.710 0.839 2.469 3.850 * 1.015 2.469 1.930 6.761 * 6.929 3.850 6.761 4.344 * */ double a[] = {7.027, 8.710, 1.015, 6.929, 8.710, 0.839, 2.469, 3.850, 1.015, 2.469, 1.930, 6.761, 6.929, 3.850, 6.761, 4.344}; /* Query optimal work size */ dsytrd_(&uplo, &n, a, &lda, d, e, tau, &qwork, &lwork, &info); if (info != 0) { printf("ERROR, info = %d\n", info); return 1; } lwork = (int)qwork; work = (double *)malloc(sizeof(double) * lwork); /* Calculate Q */ dsytrd_(&uplo, &n, a, &lda, d, e, tau, work, &lwork, &info); dorgtr_(&uplo, &n, a, &lda, tau, work, &lwork, &info); /* output */ * a * 1.000000 0.000000 0.000000 0.000000 * 0.000000 -0.779341 0.346782 -0.521890 * 0.000000 -0.090819 -0.886612 -0.453510 * 0.000000 -0.619983 -0.306042 0.722467
Fortran Interface:
CHARACTER :: uplo = "L" PARAMETER (n = 4) PARAMETER (lda = 4) INTEGER :: info = 0 REAL(8) :: a(lda, n) REAL(8) :: d(n) REAL(8) :: e(n-1) REAL(8) :: tau(n) REAL(8), ALLOCATABLE :: work(:) REAL(8) :: qwork(1) INTEGER :: lwork = -1 DATA a / 7.027, 8.710, 1.015, 6.929, & 8.710, 0.839, 2.469, 3.850, & 1.015, 2.469, 1.930, 6.761, & 6.929, 3.850, 6.761, 4.344 / DATA tau / 1.003949, 1.125229, 1.978923, 0.000000 / EXTERNAL DSYTRD,DORGTR CALL DSYTRD(uplo, n, a, lda, d, e, tau, qwork, lwork, info) IF (info.NE.0) THEN CALL EXIT(1) END IF lwork = INT(qwork(1)) ALLOCATE(work(lwork)) CALL DSYTRD(uplo, n, a, lda, d, e, tau, qwork, lwork, info) IF (info.NE.0) THEN CALL EXIT(1) END IF CALL DORGTR(uplo, n, a, lda, tau, work, lwork, info) DEALLOCATE(work); * * Output: * a * 1.000000 0.000000 0.000000 0.000000 * 0.000000 -0.779341 0.346782 -0.521890 * 0.000000 -0.090819 -0.886612 -0.453510 * 0.000000 -0.619983 -0.306042 0.722467