- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
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[0]); coefs = (double*)mxGetPr(prhs[1]); rw2 = (double)mxGetScalar(prhs[2]); /* phi1,tau */ phi1 = coefs[0]; phi1square = pow(phi1,2); tau = coefs[1]; /* make pointers to output data */ plhs[0] = mxCreateDoubleMatrix(1,dim,mxREAL); w = mxGetPr(plhs[0]); plhs[1] = mxCreateDoubleMatrix(dim,dim,mxREAL); U = mxGetPr(plhs[1]); /* Memory */ d = (double*)mxMalloc(dim*sizeof(double)); e = (double*)mxMalloc(dim*sizeof(double)); /* e, d */ e[0:dim:1] = -phi1*tau; d[0] = 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" ); } }
Link Copied
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
vl, vu
If range = 'V', the lower and upper bounds of the interval to be searched for eigenvalues. Constraint: vl<vu.

- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page