1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
|
/*
* Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
* Copyright (C) 2008 - INRIA - Arnaud TORSET
*
* This file must be used under the terms of the CeCILL.
* This source file is licensed as described in the file COPYING, which
* you should have received as part of this distribution. The terms
* are also available at
* http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
*
*/
/*
This is a transcription of the fortran subroutine cortr.f
*/
#include "logm_internal.h"
#include "conj.h"
#include "addition.h"
#include "multiplication.h"
#include "division.h"
void cortr(doubleComplex* in, int n, int low, int high, doubleComplex* ort, doubleComplex* z){
int i=0, j=0, k=0, ii=0, iend=0, ip1=0;
double norm;
doubleComplex s;
/*.......... initialize eigenvector matrix ..........*/
for (i=1;i<=n;i++){
for (j=1;j<=n;j++){
if (i==j) z[(j-1)*n+i-1]=DoubleComplex(1,0);
else z[(j-1)*n+i-1]=DoubleComplex(0,0);
}
}
/*.......... form the matrix of accumulated transformations
from the information left by corth ..........*/
iend = high -low -1;
if (iend >= 0){
for (ii=1;ii<=iend;ii++){
i=high-ii;
if ( ( !((zreals(ort[i-1])==0) && (zimags(ort[i-1])==0) ))
&& ( !((zreals(in[(i-2)*n+i-1])==0) && (zimags(in[(i-2)*n+i-1])==0) ) ) ){
/* .......... norm below is negative of h formed in corth .........*/
norm = zreals(zmuls(zconjs(in[(i-2)*n+i-1]),ort[i-1]));
if (norm != 0 ){
ip1=i+1;
for(k=ip1;k<=high;k++){
ort[k-1] = DoubleComplex( zreals(in[(i-2)*n+k-1]), zimags(in[(i-2)*n+k-1]) );
}
for (j=i;j<=high;j++){
s=DoubleComplex(0,0);
for(k=i;k<=high;k++){
s=zadds(s, zmuls(zconjs(ort[k-1]),z[(j-1)*n+k-1]));
}
s = zrdivs(s,DoubleComplex(norm,0));
for(k=i;k<=high;k++){
z[(j-1)*n+k-1]=zadds(z[(j-1)*n+k-1], zmuls(s,ort[k-1]));
}
}
}
}
}
}
}
|
