Dev - C++ - #include vector.hpp

   
Vista:

#include vector.hpp

Publicado por David davidjcg2789@gmail.com (2 intervenciones) el 23/06/2015 06:50:10
Hola a todos!

Tengo este codigo de C++ que compila haciendo uso de #include "vector.hpp". Pero necesito poder compilarlo sin el mismo. Me gustaria que alguien me ayudara a conseguir una manera de definir de otra manera los vectores que usan esta estructura (por ejemplo Vector f y Vector txyz).

Muchas gracias.

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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
#include <math.h>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <iomanip>
#include <iterator>
#include <sstream>
#include <string>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>
#include <math.h>
#include <cmath>
#include <cstdlib>
#include <fstream>
#include <iostream>
#include <iomanip>
#include<iterator>
#include <stdio.h>
#include <cctype>
using namespace std;
#include "vector.hpp"
using namespace cpl;
//#include "Vector.hpp"
 
 
double sigma = 10; // Lorenz model constants in textbook
double b = 8.0 / 3.0;
double r = 28;
int t = 400;
Vector f(Vector txyz) {
	double t = txyz[0];
	double x = txyz[1];
	double y = txyz[2];
	double z = txyz[3];
	Vector f(4);
	f[0] = 1;
    f[1] = - sigma * x + sigma * y;
    f[2] = - x * z + r * x - y;
    f[3] = x * y - b * z;
    return f;
}
 
void RK4Step( // 4th order Runge-Kutta
Vector& y, // extended solution vector
double h) // step size
{
    Vector k1 = h * f(y);
    Vector k2 = h * f(y + 0.5 * k1);
    Vector k3 = h * f(y + 0.5 * k2);
    Vector k4 = h * f(y + k3);
    y += (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
}
 
Vector txyz(4); // global variable to hold t,x,y,z
void initialize() { // initial conditions in textbook
 
    txyz[0] = 0.0;
    txyz[1] = 0.0;
    txyz[2] = 1.0;
    txyz[3] = 0.0;
}
 
double dt = 0.001; // time step for integration
 
void PrimercruceenB() {
	char B;
	char Bi;
	double w = sqrt(b * (r - 1));
	while (txyz[1] < w){ // Cruce con plano  B en x
 
        RK4Step(txyz, dt);
    }
        // Aqui ya tengo el primer valor de x cruzando al plano B
 
        if (txyz[2] >= w) { // Cruce con plano B en y
                Bi = '3';
                cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
        }
        else {
               cout << " y < w" << endl;
                while ( txyz[2] < w) {
                        RK4Step(txyz, dt);
                }
 
                if (txyz[1] >= w) {
 
                        Bi = '3';
                       cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
                }
        }
}
 
void despuesdecruceplanoB() {
        Vector txyzP = txyz; //vector presente para comprarar con el siguiente
                while (txyz[1] >= txyzP[1]) {
                        txyzP = txyz;
                        RK4Step(txyz, dt);
                }
                cout <<"x se devuelve del plano B  " << "x = " << txyz[1] << "y = " << txyz[2] << "z = " << txyz[3] << endl; // x comienza devolverse
}
 
void PrimercruceenA() {
	char Bi;
	char A;
	Vector u;
	//double g;
	//int h;
	//const int tamano_entrada = 16;
	//double delta;
	//double in[tamano_entrada] = {20.5722, 16.9518, 11.4379, 17.7182, 13.1959, 16.9518, 16.6336, 13.1959, 20.5722, 20.5722, 16.9518, 16.6336, 13.1959, 20.5722, 16.9518, 16.6336};
	double w = sqrt(b * (r - 1));
	while (txyz[1] > -w) { // Busco cruce con plano A en x
		u = txyz;
		RK4Step(txyz, dt);
	}
	/*cout << u[3] << " " << txyz[3] << endl;
				
				delta = in[h] - u[3];
           		 h = h + 1;
				u[3] = u[3] + delta;*/
	if (txyz[2] <= -w) {
		Bi = '2';
		cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
	}
	else {
		cout << " y es > w" << endl;
		while (txyz[2] > -w) {
			RK4Step(txyz, dt);
		}
		if (txyz[1] <= -w) {
			Bi = '2';
			cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
		}
	}
}
 
void despuesdecruceplanoA () {
        Vector txyzPra = txyz; //vector presente para comprarar con el siguiente
                                while (txyz[1] <= txyzPra[1]) {
                                        txyzPra = txyz;
                                        RK4Step(txyz, dt);
 
					}
					cout <<"x se devuelve del plano A  " << "x = " << txyz[1] << "y = " << txyz[2] << "z = " << txyz[3] << endl; // x comienza devolverse
}
 
 
int main(int argc, char *argv[]) {
	int n = 6;
	int i, c;
	int h = 0;
	int j = 0;
	const int tamano_entrada = 14;
	char sim[tamano_entrada];
	double delta;
	double equis[tamano_entrada + 1];
	double w = sqrt(b * (r - 1));
	double out[tamano_entrada-1];
	//Vector u = txyz;
	char Bi;
	char A;
	Vector u;
	//double in[tamano_entrada];
 
 
	if (n==6) {
		//double in[] = {20.5722, 20.5722, 20.5722, 16.9518, 16.6336, 13.1959, 20.5722, 16.9518, 16.6336, 13.1959, 16.9518, 11.4379, 17.7182, 13.1959, 16.9518, 11.4379, 20.7618, 17.7182, 13.1959};
		 double in[] = {17.1772, 12.7452, 19.7595, 18.2817, 15.8382, 14.2816, 17.6427, 14.2049, 17.6369, 14.549, 16.9111, 10.6945, 22.135, 22.135};//, 14.2816, 17.6427, 14.5687, 17.1772, 13.2075, 18.8481, 16.9111, 10.6945};
		 initialize();
 
 
		while(txyz[0] < 50) {
			RK4Step(txyz, dt);
		}
 
		while (h < tamano_entrada + n - 2) {
			PrimercruceenB();
			cout << txyz[3] << endl; // 14.4911 que me genera 322 
			despuesdecruceplanoB();
 
			// Aqui ya se hizo la perturbacion del sistema
 
			Vector txyzPr = txyz; //vector presente para comprarar con el siguiente
                while (txyz[1] <= txyzPr[1]) { // Se sale del while cuando se regresa al plano B
                        txyzPr = txyz;
                        RK4Step(txyz, dt);
 
 
                        if (txyz[1] * txyzPr[1] < 0) {
                                cout <<"x cambio de signo" << endl; // x se dirige al siguiente plano
                                 goto Parte2;
						}
				}
				if (txyz[1] > 0) {
					goto Inicio;
					break;
				}
 
			Parte2:
				{
 
			//cout << "ENTRO A PARTE 2" << endl;
 
				if (h > tamano_entrada + n - 2) {
					break;
				}
			//PrimercruceenA(); // Aqui
 
			while (txyz[1] > -w) { // Busco cruce con plano A en x
				//u = txyz;
				RK4Step(txyz, dt);
			}
			cout << txyz[3] << endl;
			delta = in[h] - txyz[3];
			cout << "delta es " << delta << endl;
			h = h + 1;
			txyz[3] = txyz[3] + delta;
			cout << "zmean es " << txyz[3] << endl;
			//RK4Step(u, dt);
			if (txyz[2] <= -w) {
				Bi = '2';
				cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
			}
            equis[j] = txyz[1];
            cout << "x es " << equis[j] << endl;
            j = j + 1;
 
			despuesdecruceplanoA ();
			Vector txyzPraa = txyz;
			while (txyz[1] >= txyzPraa[1]) {
				txyzPraa = txyz;
				RK4Step(txyz, dt);
 
                if (txyz[1] * txyzPraa[1] < 0) {
                cout <<"x cambio de signo" << endl; // x se dirige al siguiente plano
                goto Inicio;
 
				} // Cruce con el plano A
			}
			if (txyz[1] < 0) {
                goto Parte2;
			}
		}
 
		Inicio:
			{
				if (h > tamano_entrada + n - 2) {
				break;
				}
				while (txyz[1] < w){ // Cruce con plano  B en x
					RK4Step(txyz, dt);
    			}
				cout << txyz[3] << endl;
				delta = in[h] - txyz[3];
				h = h + 1;
				cout << "delta es " << delta << endl;
				txyz[3] = txyz[3] + delta;
				cout << "zmean es " << txyz[3] << endl;
				if (txyz[2] >= w) { // Cruce con plano B en y
                	Bi = '3';
                	cout << "t = " << txyz[0]<< " x = " << txyz[1] << ", y = " << txyz[2] << ", z = " << txyz[3] << ",  B = " << Bi << endl;
        		}
				equis[j] = txyz[1];
				cout << "x es " << equis[j] << endl;
				j = j + 1;
 
				despuesdecruceplanoB();
			Vector txyzPr2 = txyz; //vector presente para comprarar con el siguiente
                while (txyz[1] <= txyzPr2[1]) { // Se sale del while cuando se regresa al plano B
                        txyzPr2 = txyz;
                        RK4Step(txyz, dt);
 
 
                        if (txyz[1] * txyzPr2[1] < 0) {
                                cout <<"x cambio de signo" << endl; // x se dirige al siguiente plano
                                goto Parte2;
						}
 
				}
				if (txyz[1] > 0) {
				goto Inicio;
			}
		}
		}
		for (i=n-1; i< tamano_entrada + n - 1; i++) {
		if (equis[i] < 0) {
			sim[i] = '2';
		}
		if (equis[i] > 0) {
			sim[i] = '3';
		}
 
		cout << sim[i];
	}
	}
	cout << endl;
 
	char desim = sim[n-1];
	cout << "desim es " << desim << endl;
	for (i=n; i<tamano_entrada; i++) {
		int l = 0;
		if (sim[i] == desim) {
			desim = sim[i];
			out[l] = 1;
 
		}
		if (sim[i] != desim) {
			desim = sim[i];
			out[l] = 0;
		}
		cout << out[l];
		l = l + 1;
	}
 
}
Valora esta pregunta
Me gusta: Está pregunta es útil y esta claraNo me gusta: Está pregunta no esta clara o no es útil
0
Responder

#include vector.hpp

Publicado por Rene gar (50 intervenciones) el 23/06/2015 17:33:09
La primer cuestion es porque combinas librerias de c y c++ math.h y cmath ademas que algunos estan repetidos, puedes empezar corrigiendo eso.

Hablando de #include vector.hpp te estar refiriendo a la clase vector de la stl si es asi la deficion correcta es #include<vector> pero me da mas la impresion que es un libreria que tu te creaste y estas incluyendo su implementacion
si es la segunda opcion te invito a publicar tanto la cabezera .h como la implentacion .cpp


Te invito a publicarlo.
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar

#include vector.hpp

Publicado por David (2 intervenciones) el 23/06/2015 18:34:17
Hola Rene

Gracias por tu respuesta. La parte del codigo que yo hice en realidad, es la que empieza en la linea 68, lo anterior yo lo consegui y es donde usan el vector.hpp.

Mas abajo envio la cabecera y la implementacion (que yo no cree). Sin embargo, lo que estoy busando es una manera de hacer que este codigo compile pero sin hacer uso de esto, es decir, otra forma de declarar estos vectores sin utilizar esa libreria.

De nuevo, gracias.

Vector.hpp
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
#ifndef CPL_VECTOR_HPP
#define CPL_VECTOR_HPP
 
#include <complex>
#include <iostream>
 
namespace cpl {
 
    class Vector {
        public:
 
            Vector ( int dim = 1 );
 
            Vector ( const Vector &dv );
 
            ~Vector() { delete [] v; }
 
            int dimension() const { return dim; }
 
            void resize ( const int );
 
            const double operator[] ( const int i ) const { return v[i]; }
 
            double &operator[] ( const int i ) { return v[i]; }
 
            Vector &operator = ( const Vector &dv );
 
            Vector &operator += ( const Vector &dv );
 
            Vector &operator -= ( const Vector &dv );
 
            Vector &operator *= ( double d );
 
            Vector &operator /= ( double d );
 
            double abs();
 
            double norm();
 
            double dot ( const Vector &dv );
 
            friend std::ostream &operator<< ( std::ostream &os, const Vector &dv );
 
        private:
            int dim;
            double *v;
    };
 
    inline Vector operator + ( const Vector &dv ) {
        return dv;
    }
 
    extern Vector operator - ( const Vector &dv );
 
    extern Vector operator * ( const Vector &dv, double d );
 
    extern Vector operator * ( double d, const Vector &dv );
 
    extern Vector operator / ( const Vector &dv, double d );
 
    extern Vector operator + ( const Vector &v1, const Vector &v2 );
 
    extern Vector operator - ( const Vector &v1, const Vector &v2 );
 
    class ComplexVector {
        public:
            ComplexVector ( int dim = 1 );
            ComplexVector ( const ComplexVector &cv );
            ~ComplexVector() { delete [] v; }
 
            int dimension() const { return dim; }
            const std::complex<double> operator[] ( const int i ) const { return v[i]; }
            std::complex<double> &operator[] ( const int i ) { return v[i]; }
            ComplexVector &operator = ( const ComplexVector &cv );
 
        private:
            int dim;
            std::complex<double> *v;
    };
 
    class FFT {
        public:
            FFT() { N = 0; f = 0; inverse = false; }
            void transform ( ComplexVector &data );
            void inverseTransform ( ComplexVector &data );
            Vector power ( ComplexVector &data );
 
        private:
            int N;
            ComplexVector *f;
            bool inverse;
 
            void bitReverse();
            void DanielsonLanczos ( int n );
    };
 
    extern void solveTridiagonal (
        ComplexVector &a, ComplexVector &b, ComplexVector &c,
        ComplexVector &r, ComplexVector &u );
 
    extern void solveTridiagonalCyclic (
        ComplexVector &a, ComplexVector &b, ComplexVector &c,
        std::complex<double> alpha, std::complex<double> beta,
        ComplexVector &r, ComplexVector &x );
 
}  /* end namespace cpl */
 
#endif  /* CPL_VECTOR_HPP */



Vector.cpp
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
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
#include "Vector.hpp"
 
namespace cpl {
 
    Vector::Vector ( int dim ) {
        v = new double [this->dim = dim];
        for ( int i = 0; i < dim; i++ ) v[i] = 0;
    }
 
    Vector::Vector ( const Vector &dv ) {
        v = new double [dim = dv.dim];
        for ( int i = 0; i < dim; i++ ) v[i] = dv.v[i];
    }
 
    void Vector::resize ( const int dimension ) {
        delete [] v;
        v = new double [dim = dimension];
        for ( int i = 0; i < dim; i++ ) v[i] = 0;
    }
 
    Vector &Vector::operator = ( const Vector &dv ) {
        if ( this != &dv ) {
            if ( dim != dv.dim ) {
                delete [] v;
                v = new double [dim = dv.dim];
            }
            for ( int i = 0; i < dim; i++ ) v[i] = dv[i];
        }
        return *this;
    }
 
    Vector &Vector::operator += ( const Vector &dv ) {
        for ( int i = 0; i < dim; i++ ) v[i] += dv[i];
        return *this;
    }
 
    Vector &Vector::operator -= ( const Vector &dv ) {
        for ( int i = 0; i < dim; i++ ) v[i] -= dv[i];
        return *this;
    }
 
    Vector &Vector::operator *= ( double d ) {
        for ( int i = 0; i < dim; i++ ) v[i] *= d;
        return *this;
    }
 
    Vector &Vector::operator /= ( double d ) {
        for ( int i = 0; i < dim; i++ ) v[i] /= d;
        return *this;
    }
 
    Vector operator - ( const Vector &dv ) {
        int dim = dv.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = -dv[i];
        return temp;
    }
 
    Vector operator * ( const Vector &dv, double d ) {
        int dim = dv.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] * d;
        return temp;
    }
 
    Vector operator * ( double d, const Vector &dv ) {
        int dim = dv.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] * d;
        return temp;
    }
 
    Vector operator / ( const Vector &dv, double d ) {
        int dim = dv.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = dv[i] / d;
        return temp;
    }
 
    Vector operator + ( const Vector &v1, const Vector &v2 ) {
        int dim = v1.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = v1[i] + v2[i];
        return temp;
    }
 
    Vector operator - ( const Vector &v1, const Vector &v2 ) {
        int dim = v1.dimension();
        Vector temp ( dim );
        for ( int i = 0; i < dim; i++ ) temp[i] = v1[i] - v2[i];
        return temp;
    }
 
    double Vector::abs() {
        return std::sqrt ( norm() );
    }
 
    double Vector::norm() {
        double sum = 0;
        for ( int i = 0; i < dim; i++ ) sum += v[i] * v[i];
        return sum;
    }
 
    double Vector::dot ( const Vector &dv ) {
        double sum = 0;
        for ( int i = 0; i < dim; i++ ) sum += v[i] * dv[i];
        return sum;
    }
 
    std::ostream &operator<< ( std::ostream &os, const Vector &dv ) {
        for ( int i = 0; i < dv.dim; i++ ) {
            os << dv.v[i];
            if ( i < dv.dim - 1 )
                os << '\t';
            else
                os << '\n';
        }
        return os;
    }
 
    // ComplexVector implementation
 
    ComplexVector::ComplexVector ( int dim ) {
        v = new std::complex<double> [this->dim = dim];
        for ( int i = 0; i < dim; i++ ) v[i] = 0.0;
    }
 
    ComplexVector::ComplexVector ( const ComplexVector &cv ) {
        v = new std::complex<double> [dim = cv.dim];
        for ( int i = 0; i < dim; i++ ) v[i] = cv.v[i];
    }
 
    ComplexVector &ComplexVector::operator = ( const ComplexVector &cv ) {
        if ( this != &cv ) {
            if ( dim != cv.dim ) {
                delete [] v;
                v = new std::complex<double> [dim = cv.dim];
            }
            for ( int i = 0; i < dim; i++ ) v[i] = cv[i];
        }
        return *this;
    }
 
    // FFT implementation
 
    void FFT::transform ( ComplexVector &data ) {
        N = data.dimension();
        f = &data;
        bitReverse();
        for ( int n = 1; n < N; n *= 2 )
            DanielsonLanczos ( n );
        for ( int i = 0; i < N; ++i )
            ( *f ) [i] /= std::sqrt ( double ( N ) );
    }
 
    void FFT::inverseTransform ( ComplexVector &data ) {
        inverse = true;
        transform ( data );
        inverse = false;
    }
 
    void FFT::bitReverse() {
        int j = 1;
        for ( int i = 1; i < N; ++i ) {
            if ( i < j ) {
                std::complex<double> temp = ( *f ) [i - 1];
                ( *f ) [i - 1] = ( *f ) [j - 1];
                ( *f ) [j - 1] = temp;
            }
            int k = N / 2;
            while ( k < j ) {
                j -= k;
                k /= 2;
            }
            j += k;
        }
    }
 
    void FFT::DanielsonLanczos ( int n ) {
        const double pi = 4 * atan ( 1.0 );
        std::complex<double> W ( 0, pi / n );
        W = inverse ? std::exp ( -W ) : std::exp ( W );
        std::complex<double> W_j ( 1, 0 );
        for ( int j = 0; j < n; ++j ) {
            for ( int i = j; i < N; i += 2 * n ) {
                std::complex<double> temp = W_j * ( *f ) [n + i];
                ( *f ) [n + i] = ( *f ) [i] - temp;
                ( *f ) [i] += temp;
            }
            W_j *= W;
        }
    }
 
    Vector FFT::power ( ComplexVector &data ) {
        Vector P ( 1 + N / 2 );
        P[0] = std::norm ( data[0] ) / double ( N );
        for ( int i = 1; i < N / 2; i++ )
            P[i] = ( std::norm ( data[i] ) + std::norm ( data[N - i] ) ) / double ( N );
        P[N / 2] = std::norm ( data[N / 2] ) / double ( N );
        return P;
    }
 
    // Solving tridiagonal complex matrices
 
    void solveTridiagonal (
        ComplexVector &a, ComplexVector &b, ComplexVector &c,
        ComplexVector &r, ComplexVector &u ) {
        int n = a.dimension();
        ComplexVector gamma ( n );
        std::complex<double> beta = b[0];
        u[0] = r[0] / beta;
        for ( int j = 1; j < n; j++ ) {
            gamma[j] = c[j - 1] / beta;
            beta = b[j] - a[j] * gamma[j];
            u[j] = ( r[j] - a[j] * u[j - 1] ) / beta;
        }
        for ( int j = n - 2; j >= 0; j-- )
            u[j] -= gamma[j + 1] * u[j + 1];
    }
 
    void solveTridiagonalCyclic (
        ComplexVector &a, ComplexVector &b, ComplexVector &c,
        std::complex<double> alpha, std::complex<double> beta,
        ComplexVector &r, ComplexVector &x ) {
        int n = a.dimension();
        ComplexVector bb ( n ), u ( n ), z ( n );
        std::complex<double> gamma = -b[0];
        bb[0] = b[0] - gamma;
        bb[n - 1] = b[n - 1] - alpha * beta / gamma;
        for ( int i = 1; i < n - 1; i++ )
            bb[i] = b[i];
        solveTridiagonal ( a, bb, c, r, x );
        u[0] = gamma;
        u[n - 1] = alpha;
        for ( int i = 1; i < n - 1; i++ )
            u[i] = 0.0;
        solveTridiagonal ( a, bb, c, u, z );
        std::complex<double> fact = x[0] + beta * x[n - 1] / gamma;
        fact /= 1.0 + z[0] + beta * z[n - 1] / gamma;
        for ( int i = 0; i < n; i++ )
            x[i] -= fact * z[i];
    }
 
}  /* end namespace cpl */
Valora esta respuesta
Me gusta: Está respuesta es útil y esta claraNo me gusta: Está respuesta no esta clara o no es útil
0
Comentar