1 /**********************************************************************
2 lapack_dstedc2.c:
3
4 lapack_dstedc2.c is a subroutine to find eigenvalues and eigenvectors
5 of tridiagonlized real matrix using lapack's routines dstedc.
6
7 Log of lapack_dstedc2.c:
8
9 Dec/24/2004 Released by T.Ozaki
10
11 ***********************************************************************/
12
13 #include <stdio.h>
14 #include <math.h>
15 #include <stdlib.h>
16 #include "openmx_common.h"
17 #include "lapack_prototypes.h"
18 #include "mpi.h"
19
20
lapack_dstedc2(INTEGER N,double * D,double * E,double * W,dcomplex ** ev)21 void lapack_dstedc2(INTEGER N, double *D, double *E, double *W, dcomplex **ev)
22 {
23 int i,j;
24 char *COMPZ="I";
25 INTEGER IL,IU;
26 double ABSTOL_GR=LAPACK_ABSTOL;
27 double *Z;
28 INTEGER LDZ;
29 INTEGER LWORK;
30 double *WORK;
31 INTEGER *IWORK;
32 INTEGER LIWORK;
33 INTEGER INFO;
34
35 LDZ = N;
36
37 Z = (double*)malloc(sizeof(double)*LDZ*N);
38
39 LWORK = 1 + 4*N + N*N;
40 LIWORK = 3 + 5*N;
41 WORK = (double*)malloc(sizeof(double)*LWORK);
42 IWORK = (INTEGER*)malloc(sizeof(INTEGER)*LIWORK);
43
44 F77_NAME(dstedc,DSTEDC)( COMPZ, &N, D, E, Z, &LDZ, WORK, &LWORK, IWORK, &LIWORK, &INFO );
45
46 /* store eigenvectors */
47
48 for (i=0; i<N; i++) {
49 for (j=0; j<N; j++) {
50 ev[i+1][j+1].r = Z[i*N+j];
51 ev[i+1][j+1].i = 0.0;
52 }
53 }
54
55 /* shift ko by 1 */
56 for (i=N; i>=1; i--){
57 W[i]= D[i-1];
58 }
59
60 if (INFO>0) {
61 /*
62 printf("\n error in dstedc_, info=%d\n\n",INFO);
63 */
64 }
65 if (INFO<0) {
66 printf("info=%d in dstedc_\n",INFO);
67 MPI_Finalize();
68 exit(0);
69 }
70
71 free(Z);
72 free(WORK);
73 free(IWORK);
74 }
75