?sy(he)trd
将对称矩阵或Hermite矩阵通过相似变换成对称三对角T,即Q**H*A*Q=T。
接口定义
C Interface:
void ssytrd_(const char *uplo, const int *n, float *a, const int *lda, float *d, float *e, float *tau, float *work, const int *lwork, int *info);
void dsytrd_(const char *uplo, const int *n, double *a, const int *lda, double *d, double *e, double *tau, double *work, const int *lwork, int *info);
void chetrd_(const char *uplo, const int *n, float _Complex *a, const int *lda, float *d, float *e, float _Complex *tau, float _Complex *work, const int *lwork, int *info);
void zhetrd_(const char *uplo, const int *n, double _Complex *a, const int *lda, double *d, double *e, double _Complex *tau, double _Complex *work, const int *lwork, int *info);
Fortran Interface:
SSYTRD(UPLO,N,A,LDA,D,E,TAU,WORK,LWORK,INFO);
DSYTRD(UPLO,N,A,LDA,D,E,TAU,WORK,LWORK,INFO);
CHETRD(UPLO,N,A,LDA,D,E,TAU,WORK,LWORK,INFO);
ZHETRD(UPLO,N,A,LDA,D,E,TAU,WORK,LWORK,INFO);
参数
参数名 |
类型 |
描述 |
输入/输出 |
|---|---|---|---|
uplo |
字符型 |
|
输入 |
N |
整数型 |
矩阵A的维数,N ≥ 0。 |
输入 |
A |
|
|
输入、输出 |
lda |
整数型 |
矩阵A的主维,LDA ≥ max(1,N)。 |
输入 |
D |
|
三对角矩阵T的对角元素,即D(i)=A(i,i)。 |
输出 |
E |
|
三对角矩阵T非对角元素。
|
输出 |
tau |
|
tau(i)必须包含基本反射器的常量因子,维数为N-1。 |
输出 |
work |
|
维数(max(1,lwork))。 若info=0,work(1)返回最优lwork值。 |
输出 |
lwork |
整数型 |
数组work的维数。lwork ≥ max(1, N),对于最优性能:lwork ≥ N*NB,其中NB为最优块大小如果lwork=-1,则该例程只计算work数组的最优大小,并以work数组的第一个值返回。 |
输入 |
info |
整数型 |
|
输出 |
依赖
#include "klapack.h"
示例
C Interface:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 | 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) { return ERROR; } lwork = (int)qwork; work = (double *)malloc(sizeof(double) * lwork); /* Calculate Q */ dsytrd_(&uplo, &n, a, &lda, d, e, tau, work, &lwork, &info); free(work); /* * Output: * A output (stored in column-major) * 7.027000 -11.176107 0.051041 0.348434 * 8.710000 7.026600 6.375997 0.224161 * 1.015000 2.469000 3.358516 -1.395021 * 6.929000 3.850000 6.761000 -3.272116 * * D output (stored in column-major) * 7.027000 7.026600 3.358516 -3.272116 * * E output (stored in column-major) * -11.176107 6.375997 -1.395021 */ |
Fortran Interface:
CHARACTER :: uplo = "L"
PARAMETER (n = 4)
PARAMETER (lda = 4)
INTEGER :: info = 0
REAL(8) :: tau(4)
REAL(8) :: qwork(1)
INTEGER :: lwork = -1
REAL(8), ALLOCATABLE :: work(:)
REAL(8), ALLOCATABLE :: d(:)
REAL(8), ALLOCATABLE :: e(:)
*
* 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
*
DATA tau /1.003949, 1.125229, 1.978923, 0.000000/
REAL(8) :: a(m, n)
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 /
EXTERNAL DORGLQ
* Query optimal work size
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))
* Calculate Q
CALL DSYTRD(uplo, n, a, lda, d, e, tau, work, lwork, info)
DEALLOCATE(work)
* Output:
* A output (stored in column-major)
* 7.027000 -11.176107 0.051041 0.348434
* 8.710000 7.026600 6.375997 0.224161
* 1.015000 2.469000 3.358516 -1.395021
* 6.929000 3.850000 6.761000 -3.272116
*
* D output (stored in column-major)
* 7.027000 7.026600 3.358516 -3.272116
*
* E output (stored in column-major)
* -11.176107 6.375997 -1.395021