006
05.10.2006, 19:59 Uhr
~testo
Gast
|
hmm danke ich könnte die eine klasse mal posten in der der fehler liegen muss aber obs was bringt...?? kompilieren werdet ihr es net können weil ich fortran aufrufe drin hab....aber hier mal der teil der dafür verantwortlich sein muss:
| C++: |
std::vector< double >
Spai::Spai_Algorithm(const std::vector< double >& A,
const std::vector< double >& P,
const int dimension
)
{
int m,
n,
lda,
k,
info,
lwork,
one = 1;
char *SIDE="L",
*TRANS="T", //from left and transpose
*UPLO="U",
*NCHAR="N"; // upper triangular matrix not transposed
std::vector< double > work,
tau,
M,
A_hat,
e_k_hat,
p_col(dimension);
std::vector< double >::iterator iter;
for(int col = 0; col < dimension; col++)
{
Get_Pattern_Column (p_col,
P,
col,
dimension);
Create_e_k (col,
dimension,
e_k_hat);
Create_A_Hat__e_k_Hat (A,
n,
m,
p_col,
A_hat,
e_k_hat);
Update_Lapack_Values (m,
n,
k,
lda,
tau,
lwork,
work);
get_opt_lwork (A_hat,
m,
n);
dgeqrf_ (&m,
&n,
&A_hat[0],
&lda,
&tau[0],
&work[0],
&lwork,
&info);
//Computing : Q^T * e_k_hat
dormqr_ (SIDE,
TRANS,
&m,
&n,
&k,
&A_hat[0],
&lda,
&tau[0],
&e_k_hat[0],
&lda,
&work[0],
&lwork,
&info);
Ectract_quadratic_content (A_hat,
m,
n);
dtrtrs_ (UPLO,
NCHAR,
NCHAR,
&n,
&one,
&A_hat[0],
&n,
&e_k_hat[0],
&n,
&info);
Insert_m_k_hat (M,
p_col,
e_k_hat);
}
return M;
}
void
Spai::Get_Pattern_Column(std::vector< double >& p_col,
const std::vector< double >& P,
const int col,
const int dimension)
{
p_col.clear();
p_col.insert(p_col.begin(),
P.begin() + dimension * col,
P.begin() + dimension * ( col + 1 ));
}
void
Spai::Create_e_k(const int col,
const int dimension,
std::vector< double >& e_k_hat)
{
e_k_hat.clear();
e_k_hat.resize(dimension, 0.0);
e_k_hat[col] = 1.0;
}
void
Spai::Create_A_Hat__e_k_Hat(const std::vector< double >& A,
int& A_hat_dim_x,
int& A_hat_dim_y,
const std::vector< double >& p_col,
std::vector< double >& A_hat,
std::vector< double >& e_k_hat)
{
A_hat_dim_x = 0;
A_hat_dim_y = p_col.size();
A_hat.clear();
for(int y = 0; y < p_col.size(); y++)
{
if(fabs(p_col[y] - 1.0) <= epsilon)
{
A_hat.insert(A_hat.end(),
A.begin() + p_col.size() * y,
A.begin() + p_col.size() * (y + 1));
A_hat_dim_x++;
}
}
Erase_Zero_Rows(A_hat, A_hat_dim_x, A_hat_dim_y, e_k_hat);
}
void
Spai::Erase_Zero_Rows(std::vector< double >& matrix,
int& x_dim,
int& y_dim,
std::vector< double >& vector)
{
std::vector< bool > mask( y_dim );
for ( int y = 0; y < y_dim; ++y )
{
mask[ y ] = true;
if(fabs(vector[ y ] - 1.0) <= epsilon)
mask[ y ] = false;
for ( int x = 0; x < x_dim; ++x )
{
if ( matrix[ x * y_dim + y ] != 0 )
{
mask[ y ] = false;
break;
}
}
}
matrix.erase( std::remove_if( matrix.begin(),
matrix.end(),
IS_MASKED( matrix, mask ) ),
matrix.end() );
vector.erase( std::remove_if( vector.begin(),
vector.end(),
IS_MASKED( vector, mask ) ),
vector.end() );
y_dim = matrix.size() / x_dim;
}
void
Spai::Update_Lapack_Values(const int m,
const int n,
int& k,
int& lda,
std::vector< double >& tau,
int& lwork,
std::vector< double >& work)
{
tau.clear();
work.clear();
k = std::min(m,n);
lda = std::max(m,1);
lwork = std::max(n,1);
tau.resize(k);
work.resize(std::max(1,lwork));
}
void
Spai::Insert_m_k_hat(std::vector< double >& M,
const std::vector< double >& p_col,
const std::vector< double >& e_k_hat)
{
std::vector< double >::const_iterator el = e_k_hat.begin();
for( std::vector< double >::const_iterator i = p_col.begin();
i != p_col.end();
++i )
{
if( !(fabs(*i - 0.0) <= epsilon) )
M.push_back(*el++);
else
M.push_back(0.0);
}
}
void
Spai::Ectract_quadratic_content(std::vector< double >& A_hat,
const int m,
const int n)
{
if (m <= n)
A_hat.erase(A_hat.begin() + m * m, A_hat.end());
else
{
std::vector< bool > mask( m );
for ( int y = 0; y < m; ++y )
{
mask[ y ] = false;
if(y >= n)
mask[ y ] = true;
}
A_hat.erase( std::remove_if( A_hat.begin(),
A_hat.end(),
IS_MASKED( A_hat, mask ) ),
A_hat.end() );
}
}
int
Spai::get_opt_lwork(std::vector< double >& A_hat,
int m,
int n)
{
int lw = -1,
info = 0,
lwork = 0;
double work = 0.0;
std::vector< double > tau;
dgeqrf_(&m,
&n,
&A_hat[0],
&m,
&tau[0],
&work,
&lw,
&info);
lwork = static_cast< int > (work);
return lwork;
}
|
|
