/*--------------------------------------------------------------------------*/
/*--------------------------------------------------------------------------*/
-// dgemm performs one of the matrix-matrix operations
-extern int C2F(dgemm)(char* transa, char* transb, int* m, int* n, int* k,
- double* alpha, double* A, int* lda, double* B, int* ldb,
- double* beta, double* C, int* ldc);
-// zgemm performs one of the matrix-matrix operations
-extern int C2F(zgemm)(char* transa, char* transb, int* m, int* n, int* k,
- doublecomplex* alpha, doublecomplex* A, int* lda, doublecomplex* B,
- int* ldb, doublecomplex* beta, doublecomplex* C, int* ldc);
-
// dgemv performs matrix-vector operations
extern int C2F(dgemv) (char* trans, int* m, int* n, double* alpha, double* A, int* lda,
double* x, int* incx, double* beta, double* y, int* incy);
int * lworkl, double * rwork, int * info);
/*--------------------------------------------------------------------------*/
+/*--------------------------------------------------------------------------*/
+extern int C2F(dsymv)(char* UPLO, int* N, double* ALPHA, double* A, int* LDA, double* X, int* INCX, double* BETA, double* Y, int* INCY);
+/*--------------------------------------------------------------------------*/
+
+/*--------------------------------------------------------------------------*/
+extern int C2F(daxpy)(int* N, double* DA, double* DX, int* INCX, double* DY, int* INCY);
+/*--------------------------------------------------------------------------*/
+
+/*--------------------------------------------------------------------------*/
+extern int C2F(zaxpy)(int* N, doublecomplex* ZA, doublecomplex* ZX, int* INCX, doublecomplex* ZY, int* INCY);
+/*--------------------------------------------------------------------------*/
+
static double alpha = 1.;
static double beta = 0.;
double* WORKEV = NULL;
doublecomplex* WORKEVC = NULL;
- doublecomplex mSIGMA = {.r = -SIGMA.r, .i = -SIGMA.i };
+ doublecomplex mSIGMA = { -SIGMA.r, -SIGMA.i };
double* R = NULL;
doublecomplex* RC = NULL;
if (DI[i] == 0)
{
C2F(dgemv) ("n", &N, &N, &alpha, AR, &N, Z + N * i, &iOne, &beta, temp1, &iOne);
- eigenvalue[i] = (doublecomplex)
- {
- C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne), 0
- };
+ eigenvalue[i].r = C2F(ddot) (&N, Z + N * i, &iOne, temp1, &iOne);
+ eigenvalue[i].i = 0;
i = i + 1;
}
else
C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp2, &iOne);
imag_part = C2F(ddot) (&N, Z + N * i, &iOne, temp2, &iOne) - \
C2F(ddot) (&N, Z + N * (i + 1), &iOne, temp1, &iOne);
- eigenvalue[i] = (doublecomplex)
- {
- real_part, imag_part
- };
- eigenvalue[i + 1] = (doublecomplex)
- {
- real_part, -imag_part
- };
+ eigenvalue[i].r = real_part;
+ eigenvalue[i].i = imag_part;
+ eigenvalue[i + 1].r = real_part;
+ eigenvalue[i + 1].i = -imag_part;
i = i + 2;
}
}
{
for (i = 0; i < nev + 1; i++)
{
- eigenvalue[i] = (doublecomplex)
- {
- DR[i], DI[i]
- };
+ eigenvalue[i].r = DR[i];
+ eigenvalue[i].i = DI[i];
}
}
{
for (j = 0; j < N; j++)
{
- eigenvector[i * N + j] = (doublecomplex)
- {
- Z[i * N + j], Z[(i + 1) * N + j]
- };
- eigenvector[(i + 1) * N + j] = (doublecomplex)
- {
- Z[i * N + j], -Z[(i + 1) * N + j]
- };
+ eigenvector[i * N + j].r = Z[i * N + j];
+ eigenvector[i * N + j].i = Z[(i + 1) * N + j];
+ eigenvector[(i + 1) * N + j].r = Z[i * N + j];
+ eigenvector[(i + 1) * N + j].i = -Z[(i + 1) * N + j];
}
+
i = i + 2;
}
else
{
for (j = 0; j < N; j++)
{
- eigenvector[i * N + j] = (doublecomplex)
- {
- Z[i * N + j], 0
- };
+ eigenvector[i * N + j].r = Z[i * N + j];
+ eigenvector[i * N + j].i = 0;
}
+
i = i + 1;
}
}