Community
cancel
Showing results for 
Search instead for 
Did you mean: 
A___Fiori
Beginner
29 Views

problem with LAPACKE_dstemr

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" );
    }
}

 

0 Kudos
1 Reply
Gennady_F_Intel
Moderator
29 Views

vl, vu

If range = 'V', the lower and upper bounds of the interval to be searched for eigenvalues. Constraint: vl<vu.