Mercurial > matrix-functions

comparison funm_files/swap.c @ 0:8f23314345f4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Create local repository for matrix toolboxes. Step #0 done.
author Antonio Pino Robles <data.script93@gmail.com>
date Wed, 06 May 2015 14:56:53 +0200
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:8f23314345f4
1 /*
2 * C mex file for MATLAB that implements LAPACK ctrexc_ for
3 * reordering the complex Schur decomposition A = UTU^* such that
4 * T is in block form. This program combines Cswap and Crealswap
5 * which deal with the complex and the real case respectively.
6 *
7 * Input matrices 'U' and 'T' are those from Schur decompositon
8 *
9 * Called by [U,T] = swap(U,T,M)
10 *
11 * The matrix M is produced using the m-files
12 * >> m = blocking(T,noplot,delta); where delta is some tolerance
13 * default: delta = 0.1;
14 * >> [n,M,ind,totalswaps] = swapping5(m);
15 *
16 * Output matrices 'U' and 'T' are the updated Schur decomposition
17 *
18 */
19
20 #include "mex.h"
21 #include "matrix.h"
22 #include "fort.h" /* defines mat2fort and fort2mat */
23
24 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
25 {
26
27 /* compq=V then the matrix of Schur vectors is updated */
28
29 char *compq = "V", msg[80];
30 int n, ldt, ldq, i, info, ifst, ilst, lengthm;
31 double *t, *q, *m, *work;
32 mxArray *mm;
33
34 /* expect 3 inputs and 2 outputs */
35 if ((nrhs != 3) || (nlhs != 2)){
36 mexErrMsgTxt("Expected 3 inputs and 2 outputs");
37 }
38
39 /* Sizes of stuff */
40
41 n = mxGetN( prhs[0] );
42 ldt = n;
43 ldq = n;
44 lengthm = mxGetM( prhs[2] );
45
46 /* Copy input and output matrices to stuff */
47
48 mm = mxDuplicateArray( prhs[2] );
49 m = mxGetPr(mm);
50
51 if (mxIsComplex( prhs[1] )) {
52
53 /* Convert complex input to complex FORTRAN format */
54
55 q = mat2fort( prhs[0], ldq, n );
56 t = mat2fort( prhs[1], ldt, n ); }
57
58 else {
59
60 plhs[0] = mxCreateDoubleMatrix(n,n,mxREAL);
61 plhs[1] = mxCreateDoubleMatrix(n,n,mxREAL);
62
63 plhs[1] = mxDuplicateArray( prhs[1] );
64 t = mxGetPr(plhs[1]);
65 plhs[0] = mxDuplicateArray( prhs[0] );
66 q = mxGetPr(plhs[0]);
67 }
68
69 /* Allocate workspace */
70
71 work = (double *)mxCalloc(n,sizeof(double));
72
73 /* Do the loop */
74
75 for ( i = 0; i < lengthm; i++) {
76 info = 0;
77 ifst = m[lengthm + i];
78 ilst = m[i];
79 if (mxIsComplex( prhs[1] )) {
80 ztrexc(compq,&n,t,&ldt,q,&ldq,&ifst,&ilst,&info); }
81 else {
82 dtrexc(compq,&n,t,&ldt,q,&ldq,&ifst,&ilst,work,&info);
83 }
84 if (info < 0){
85 sprintf(msg, "The routine DTREXC has detected an error");
86 mexErrMsgTxt(msg);
87 }
88 }
89
90 /* Convert output to MATLAB format */
91
92 if (mxIsComplex( prhs[1] )) {
93 plhs[0] = fort2mat( q, ldq, ldq, n );
94 plhs[1] = fort2mat( t, ldt, ldt, n );
95 }
96
97 /* Free up memory */
98
99 mxFree(work);
100 }
101