Intel® oneAPI Math Kernel Library
Ask questions and share information with other developers who use Intel® Math Kernel Library.

## problem with LAPACKE_dstemr Beginner
142 Views

Hello! I am trying to use the function 'LAPACKE_dstemr' with a symmetric tridiagonal matrix 'T' but my code crashes.  The  matrix T (dimxdim) in general has the following form:

| 1.0        a         0         0         0   |
| a           b         a         0         0   |
T =   | 0           a         b         a         0   | *tau
| 0           0         a         b          a  |
| 0           0         0         a         1.0|

where a = -phi1, b=(1.0 + phi1^2)

Let  for example that coefs=[phi1 tau]=[0.7 1] and  dim=200. I want to compute all the eigenvalues and eigenvectors. I cant find an example with this function, please could you tell me what I am doing wrong. My mex-code is the following. Thank you very much.

```#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <string.h> /* is needed for the function 'strcmp' */

/* Auxiliary routines prototypes */
extern void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda );
extern double sprod(int dim, double *a);

/* main fucntion */
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{
int i,dim,info,n;
int tryrac = true;
double *coefs,phi1,tau,phi1square,rw2,*U,*w,*e,*d;

/* make pointers to input data */
dim = (int)mxGetScalar(prhs);
coefs = (double*)mxGetPr(prhs);
rw2 = (double)mxGetScalar(prhs);

/* phi1,tau */
phi1 = coefs;
phi1square = pow(phi1,2);
tau = coefs;

/* make pointers to output data */
plhs = mxCreateDoubleMatrix(1,dim,mxREAL);
w = mxGetPr(plhs);

plhs = mxCreateDoubleMatrix(dim,dim,mxREAL);
U = mxGetPr(plhs);

/* Memory */
d = (double*)mxMalloc(dim*sizeof(double));
e = (double*)mxMalloc(dim*sizeof(double));

/* e, d */
e[0:dim:1] = -phi1*tau;
d = tau;
d[dim-1] = tau;
d[1:dim-2:1] = (1.0 + phi1square)*tau;

print_matrix( "e", 1, dim, e, 1 );
print_matrix( "d", 1, dim, d, 1 );

/* Solve eigenproblem */
n = dim;
info = LAPACKE_dstemr(LAPACK_COL_MAJOR, 'V', 'A', dim,
d, e, 0, 0, 0, 0, &n,  w, U, dim, dim, &n, &tryrac );

mexPrintf( "n=%d\t dim=%d\n", n,dim );

/* release memory */
mxFree(e),mxFree(d);
return ;
}

/* Auxiliary routine: printing a matrix */
void print_matrix( char* desc, MKL_INT m, MKL_INT n, double* a, MKL_INT lda ) {
MKL_INT i, j;
mexPrintf( "\n %s\n", desc );
for( i = 0; i < m; i++ ) {
for( j = 0; j < n; j++ ) mexPrintf( " %6.2f", a[i+j*lda] );
mexPrintf( "\n" );
}
}
```  