You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
312 lines
6.6 KiB
312 lines
6.6 KiB
#include <cmath>
|
|
#include <iostream>
|
|
#include "array.h"
|
|
#include <fstream>
|
|
|
|
using namespace std;
|
|
|
|
|
|
template <typename T>
|
|
ArrayFP<T>::ArrayFP(unsigned int dim)
|
|
{
|
|
N = dim;
|
|
if (N == 0) entry = 0;
|
|
else
|
|
entry = new T[N];
|
|
|
|
reset();
|
|
}
|
|
|
|
template <typename T>
|
|
ArrayFP<T>::~ArrayFP(){
|
|
|
|
if (entry != 0)
|
|
delete [] entry;
|
|
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::setup(int dim)
|
|
{
|
|
N = dim;
|
|
if (N == 0) entry = 0;
|
|
else
|
|
entry = new T[N];
|
|
|
|
reset();
|
|
}
|
|
|
|
template <typename T>
|
|
T ArrayFP<T>::magnitude()
|
|
{
|
|
int i;
|
|
T mag;
|
|
|
|
mag = 0.0;
|
|
for (i = 0; i < N; i ++) {
|
|
mag += (entry[i] * entry[i]);
|
|
}
|
|
|
|
mag = sqrt(mag);
|
|
return mag;
|
|
}
|
|
|
|
template <typename T>
|
|
T ArrayFP<T>::getentry(int m) const { return entry[m]; }
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::addArray(const ArrayFP &xvtr)
|
|
{
|
|
int i;
|
|
T cval;
|
|
|
|
if (xvtr.N != N) return;
|
|
for (i = 0; i < N; i ++) {
|
|
cval = xvtr.getentry(i);
|
|
|
|
entry[i] = entry[i] + cval;
|
|
}
|
|
}
|
|
|
|
template <typename T>
|
|
T ArrayFP<T>::operator * (const ArrayFP &operand2) const
|
|
{
|
|
T cval = 0.0;
|
|
int i;
|
|
|
|
for (i = 0; i < N; i ++) {
|
|
cval = cval + entry[i] * operand2.entry[i];
|
|
}
|
|
|
|
return cval;
|
|
}
|
|
|
|
template <typename T>
|
|
ArrayFP<T>& ArrayFP<T>::operator=(const ArrayFP &right)
|
|
{
|
|
// int i;
|
|
|
|
// N = right.N;
|
|
// entry = new T[N];
|
|
// for (i = 0; i < N; i ++) entry[i] = right.entry[i];
|
|
return *this;
|
|
}
|
|
|
|
|
|
template <>
|
|
ArrayFP<float>& ArrayFP<float>::operator=(const ArrayFP<float> &right)
|
|
{
|
|
cblas_scopy(right.N, right.entry, 1, entry, 1);
|
|
return *this;
|
|
}
|
|
|
|
template <>
|
|
ArrayFP<double>& ArrayFP<double>::operator=(const ArrayFP<double> &right)
|
|
{
|
|
cblas_dcopy(right.N, right.entry, 1, entry, 1);
|
|
return *this;
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::print()
|
|
{
|
|
cout.setf(ios::scientific);
|
|
cout.precision(15);
|
|
cout << "begin"<< endl;
|
|
for (int i = 0 ; i < N ; i++)
|
|
cout << "x(" << i+1 << ") = " << entry[i] << endl;
|
|
|
|
cout << "end"<< endl;
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::print(char* name)
|
|
{
|
|
cout << name << endl;
|
|
cout.setf(ios::scientific);
|
|
cout.precision(7);
|
|
cout << "begin"<< endl;
|
|
int M=N/12;
|
|
for (int j = 0 ; j < M ; j++){
|
|
for (int i = 0 ; i < 12 ; i++){
|
|
// cout << "x(" << i+1 << ") = " << entry[j*12+i] <<" ";
|
|
cout << entry[j*12+i] << " ";
|
|
}cout << endl;
|
|
}
|
|
cout << "end"<< endl;
|
|
}
|
|
|
|
// y = alpha*x
|
|
template <typename T>
|
|
void ArrayFP<T>::ax(T alpha, ArrayFP& x)
|
|
{ }
|
|
|
|
template <>
|
|
void ArrayFP<float>::ax(float alpha, ArrayFP<float>& x)
|
|
{
|
|
cblas_scopy(N, x.entry, 1, entry, 1);
|
|
cblas_sscal(N, alpha, entry, 1);
|
|
}
|
|
|
|
template <>
|
|
void ArrayFP<double>::ax(double alpha, ArrayFP<double>& x)
|
|
{
|
|
cblas_dcopy(N, x.entry, 1, entry, 1);
|
|
cblas_dscal(N, alpha, entry, 1);
|
|
}
|
|
|
|
// y = alpha*x + y
|
|
template <typename T>
|
|
void ArrayFP<T>::axpy(T alpha, ArrayFP& x)
|
|
{
|
|
|
|
}
|
|
|
|
template <>
|
|
void ArrayFP<float>::axpy(float alpha, ArrayFP<float>& x)
|
|
{
|
|
cblas_saxpy(N, alpha, x.entry, 1, entry, 1);
|
|
}
|
|
|
|
template <>
|
|
void ArrayFP<double>::axpy(double alpha, ArrayFP<double>& x)
|
|
{
|
|
cblas_daxpy(N, alpha, x.entry, 1, entry, 1);
|
|
}
|
|
|
|
// y = alpha*x + beta*y
|
|
template <typename T>
|
|
void ArrayFP<T>::axpby(T alpha, T beta, ArrayFP& x)
|
|
{
|
|
|
|
}
|
|
|
|
template <>
|
|
void ArrayFP<float>::axpby(float alpha, float beta, ArrayFP<float>& x)
|
|
{
|
|
cblas_sscal(N, beta, entry, 1);
|
|
cblas_saxpy(N, alpha, x.entry, 1, entry, 1);
|
|
}
|
|
|
|
template <>
|
|
void ArrayFP<double>::axpby(double alpha, double beta, ArrayFP<double>& x)
|
|
{
|
|
cblas_dscal(N, beta, entry, 1);
|
|
cblas_daxpy(N, alpha, x.entry, 1, entry, 1);
|
|
}
|
|
|
|
// y = x + beta*y
|
|
template <typename T>
|
|
void ArrayFP<T>::xpby(T beta, ArrayFP& x)
|
|
{
|
|
// #pragma omp parallel for schedule(static)
|
|
// for (int i = 0; i < dim_; i++) p_[i] = x.p_[i] + beta*p_[i];
|
|
axpby(1., beta, x);
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::WriteToASCII(char* FileName, int& timeStep)
|
|
{
|
|
char Field_TimeLog[180];
|
|
// sprintf(Field_TimeLog, "%s_%d.ico", FileName, timeStep);
|
|
sprintf(Field_TimeLog, "%s", FileName);
|
|
ofstream Outfile(Field_TimeLog, ios_base::out);
|
|
Outfile.setf(ios::scientific, ios::floatfield);
|
|
Outfile.precision(15);
|
|
if(!Outfile)
|
|
cout << "ArrayFP:: Error in opening file: " << Field_TimeLog << " for write " << endl;
|
|
|
|
Outfile << timeStep << '\n';
|
|
for (int i = 0 ; i < N ; i++)
|
|
Outfile << entry[i] << '\n';
|
|
|
|
Outfile.close();
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::ReadFromASCII(char* FileName, int& timeStep)
|
|
{
|
|
char Field_TimeLog[180];
|
|
T CurrentValue = 0.0;
|
|
// sprintf(Field_TimeLog, "%s_%d.ico", FileName, timeStep);
|
|
sprintf(Field_TimeLog, "%s", FileName);
|
|
ifstream InputFile(Field_TimeLog, ios_base::in);
|
|
if (!InputFile)
|
|
cout << "ArrayFP:: Can not open file: " << Field_TimeLog << "for read" << endl;
|
|
|
|
InputFile >> timeStep;
|
|
for (int i = 0 ; i < N ; i++){
|
|
InputFile >> CurrentValue;
|
|
entry[i] = CurrentValue;
|
|
}
|
|
InputFile.close();
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::WriteToASCII(char* FileName)
|
|
{
|
|
char Field_TimeLog[180];
|
|
sprintf(Field_TimeLog, "%s", FileName);
|
|
ofstream Outfile(Field_TimeLog, ios_base::out);
|
|
Outfile.setf(ios::scientific, ios::floatfield);
|
|
Outfile.precision(15);
|
|
if(!Outfile)
|
|
cout << "ArrayFP:: Error in opening file: " << Field_TimeLog << " for write " << endl;
|
|
|
|
for (int i = 0 ; i < N ; i++)
|
|
Outfile << entry[i] << '\n';
|
|
|
|
Outfile.close();
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::ReadFromASCII(char* FileName)
|
|
{
|
|
char Field_TimeLog[180];
|
|
T CurrentValue = 0.0;
|
|
sprintf(Field_TimeLog, "%s", FileName);
|
|
ifstream InputFile(Field_TimeLog, ios_base::in);
|
|
if (!InputFile)
|
|
cout << "ArrayFP:: Can not open file: " << Field_TimeLog << "for read" << endl;
|
|
|
|
for (int i = 0 ; i < N ; i++){
|
|
InputFile >> CurrentValue;
|
|
setentry(i, CurrentValue);
|
|
}
|
|
InputFile.close();
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::WriteToBinary(char* FileName)
|
|
{
|
|
// char Field_TimeLog[180];
|
|
// sprintf(Field_TimeLog, "%s", FileName);
|
|
// ofstream Outfile(Field_TimeLog, ios_base::out, ios::binary);
|
|
// if(!Outfile)
|
|
// cout << "ArrayFP:: Error in opening file: " << Field_TimeLog << " for write " << endl;
|
|
//
|
|
// buffer = new char[size];
|
|
// T value;
|
|
// for (int i = 0 ; i < N ; i++){
|
|
// value = entry[i];
|
|
// Outfile.write(buffer, value);
|
|
// }
|
|
// Outfile.close();
|
|
}
|
|
|
|
template <typename T>
|
|
void ArrayFP<T>::ReadFromBinary(char* FileName)
|
|
{
|
|
// char Field_TimeLog[180];
|
|
// T CurrentValue = 0.0;
|
|
// sprintf(Field_TimeLog, "%s", FileName);
|
|
// ifstream InputFile(Field_TimeLog, ios_base::in, ios::binary);
|
|
// if (!InputFile)
|
|
// cout << "ArrayFP:: Can not open file: " << Field_TimeLog << "for read" << endl;
|
|
//
|
|
// for (int i = 0 ; i < N ; i++){
|
|
// InputFile.read((char *)(&CurrentValue), sizeof(CurrentValue));
|
|
// setentry(i, CurrentValue);
|
|
// }
|
|
// InputFile.close();
|
|
}
|
|
|