C++ QMat::toData方法代码示例

void RMat::QMat::SVD(QMat & U, QMat & D, QMat & V)
{
#ifdef COMPILE_IPP
// 	double ippA[rows*cols];
	double ippU[rows*rows], ippD[rows*cols], ippV[cols*cols];
	int i;
	
	U = QMat(rows, rows);
	D = QMat(rows, cols);
	V = QMat(cols, cols);

// 	for(i=0; i<rows*cols;i++)
// 	  ippA[i] = toData()[i];
	
	
	for(i=0; i<rows*rows;i++)
	  U.toData()[i] = ippU[i];
	
	for(i=0; i<cols;i+=1)
	  D.toData()[i*(cols+1)] = ippD[i];
	
	for(i=0; i<cols*cols;i++)
	  V.toData()[i] = ippV[i];
#else
	qFatal("QMat does not support eigen svd");
#endif

}

/**
 * \brief Copy function
 *
 * Return a copy of this matrix
 * @return QMat matrix copied
*/
QMat QMat::copy()
{
	QMat A (rows, cols);
#ifdef COMPILE_IPP
	ippmCopy_ma_32f_SS ( toData(), rows*cols*sizeof ( T ), cols*sizeof ( T ), sizeof ( T ), A.toData(),rows*cols*sizeof ( T ), cols*sizeof ( T ), sizeof ( T ),  cols, rows, 1 );
#else
	memcpy(A.toData(), toData(), rows*cols*sizeof(T));
#endif
	return A;
}

/**
 * \brief Change the order of the indexes in the matrix \f$ C = this^t \f$
 * @return QMat new matrix with indexes reordered
 */
QMat QMat::transpose( ) const
{
	QMat C ( cols, rows );
#ifdef COMPILE_IPP
	ippmTranspose_m_32f ( toDataConst(), cols*sizeof ( T ), sizeof ( T ), cols, rows, C.toData(), rows*sizeof ( T ), sizeof ( T ) );
#else
	for (int fila=0; fila<rows; fila++ )
		for ( int columna=0; columna<cols; columna++ )
			C(columna,fila) = operator()(fila,columna);
#endif
	return C;
}

QMat RMat::QMat::eigenValsVectors ( QVec & vals )
{
	if ( isSquare() and ( vals.size() ==cols ) )
	{
		QMat R ( rows,cols,0.);
#ifdef COMPILE_IPP
			T valsIm[cols];
			int SizeBytes;
			ippmEigenValuesVectorsGetBufSize_32f(cols, &SizeBytes);
			Ipp8u pBuffer[SizeBytes];
			IppStatus status = ippmEigenValuesVectorsRight_m_32f (toData(), cols*sizeof(T), sizeof(T), R.toData(), cols*sizeof(T), sizeof(T), vals.getWriteData(),valsIm, cols, pBuffer); 

			if ( status == ippStsNoErr )
				return R;
			else
			{
				printf ( "EigenValsVectors error: %d, %s\n", status, ippGetStatusString ( status ) );
				( *this ).print ( "this" );
				QString s="QMat::eigenVectors() - Error returned by IPP"; throw s;
			}
#else
// 			qFatal("eigen values using Eigen has not been properly tested yet");
			double data[this->getDataSize()];
			const T *dataBuffer = this->getReadData();
			for(int i = 0; i < this->getDataSize(); i++)
				data[i] = (double)dataBuffer[i];

			gsl_matrix_view m = gsl_matrix_view_array (data, rows,rows);
			gsl_vector *eval = gsl_vector_alloc (rows);
			gsl_matrix *evec = gsl_matrix_alloc (rows, rows);

			gsl_eigen_symmv_workspace * w = gsl_eigen_symmv_alloc (rows);
			gsl_eigen_symmv (&m.matrix, eval, evec, w);
			gsl_eigen_symmv_free (w);
			gsl_eigen_symmv_sort (eval, evec, GSL_EIGEN_SORT_ABS_DESC);

			for(int i=0;i<rows;i++)
				vals(i) = gsl_vector_get(eval,i);
			QMat P(evec);
			
			gsl_vector_free (eval);
			gsl_matrix_free (evec);

			return P;
#endif
	}
	QString s="QMat::EigenValsVectors() - Not Square Matrix"; throw s;
	return QMat();
}

QMat QMat::invert( ) const
{
	Q_ASSERT( isSquare() );
	QMat R ( rows, cols );
#ifdef COMPILE_IPP
	T pBuffer[cols*cols+cols];
	IppStatus status = ippmInvert_m_32f ( getReadData(), cols*sizeof ( T ), sizeof ( T ), pBuffer, R.toData(), cols*sizeof ( T ), sizeof ( T ), cols );
#else
	int status =1;
	if(determinant()!=0)
	{
		status=0;
		gsl_matrix *u = *this;
		gsl_vector *s = gsl_vector_alloc (cols);
		gsl_matrix *v = gsl_matrix_alloc (cols,cols);
		gsl_vector *work = gsl_vector_alloc (cols);
		
		gsl_matrix *X = gsl_matrix_alloc (cols,cols);
		gsl_linalg_SV_decomp_mod(u, X, v, s, work);
		gsl_matrix_free(X);

		// pasamos de vuelta:
		QMat U = u;
		QMat V = v;
		gsl_matrix_free(u);
		gsl_matrix_free(v);
		gsl_vector_free(work);
		
		QMat S = QMat::zeros(rows,cols);
		for (int i = 0; i < cols; i++)
			S(i,i) = 1 / gsl_vector_get(s, i);

		gsl_vector_free(s);
		return V * S * U.transpose();
	}
#endif
	if ( status == 0 )
		return R;
	else
	{
		QString ex= "QMat::invert() - Singular matrix";
		throw ex;
	}
}

/**
 * \brief Matrix Subtraction operator: \f$ C = this - A \f$
 * @param A Matrix subtrahend
 * @return QMat Operation result
 */
QMat QMat::operator- ( const QMat & A ) const
{
	Q_ASSERT ( equalSize ( *this, A ));
	QMat C ( rows,cols );
#ifdef COMPILE_IPP
	ippmSub_mm_32f ( toDataConst(), cols*sizeof ( T ), sizeof ( T ), A.toDataConst(), cols*sizeof ( T ), sizeof ( T ), C.toData(), cols*sizeof ( T ), sizeof ( T ), cols, rows );
#else
 	for(int r=0; r<A.rows; r++)
 		for(int c=0; c<A.cols; c++)
 			C(r, c) = A(r,c) - operator()(r, c);
#endif	
	return C;
}