00001
00005 #include <bpm/bpm_messages.h>
00006 #include <bpm/bpm_nr.h>
00007
00008 double gsl_blas_dnrm2 (const gsl_vector * X)
00009 {
00029 return cblas_dnrm2 ((int)(X->size), X->data, (int)(X->stride));
00030 }
00031
00032 double cblas_dnrm2(const int N, const double *X, const int incX) {
00033
00034
00035 double scale = 0.0;
00036 double ssq = 1.0;
00037 int i;
00038 int ix = 0;
00039
00040 if (N <= 0 || incX <= 0) {
00041 return 0;
00042 } else if (N == 1) {
00043 return fabs(X[0]);
00044 }
00045
00046 for (i = 0; i < N; i++) {
00047 const double x = X[ix];
00048
00049 if (x != 0.0) {
00050 const double ax = fabs(x);
00051
00052 if (scale < ax) {
00053 ssq = 1.0 + ssq * (scale / ax) * (scale / ax);
00054 scale = ax;
00055 } else {
00056 ssq += (ax / scale) * (ax / scale);
00057 }
00058 }
00059
00060 ix += incX;
00061 }
00062
00063 return scale * sqrt(ssq);
00064 }
00065
00066 void gsl_blas_dscal (double alpha, gsl_vector * X)
00067 {
00068 cblas_dscal ((int) (X->size), alpha, X->data, (int) (X->stride));
00069 }
00070
00071 void cblas_dscal (const int N, const double alpha, double *X, const int incX)
00072 {
00073 int i;
00074 int ix;
00075
00076 if (incX <= 0) {
00077 return;
00078 }
00079
00080 ix = OFFSET(N, incX);
00081
00082 for (i = 0; i < N; i++) {
00083 X[ix] *= alpha;
00084 ix += incX;
00085 }
00086 }
00087
00088 int gsl_blas_dgemv (CBLAS_TRANSPOSE_t TransA, double alpha, const gsl_matrix * A,
00089 const gsl_vector * X, double beta, gsl_vector * Y)
00090 {
00091 const size_t M = A->size1;
00092 const size_t N = A->size2;
00093
00094 if ((TransA == CblasNoTrans && N == X->size && M == Y->size)
00095 || (TransA == CblasTrans && M == X->size && N == Y->size))
00096 {
00097 cblas_dgemv (CblasRowMajor, TransA, (int) M, (int) N, alpha, A->data,
00098 (int) (A->tda), X->data, (int) (X->stride), beta, Y->data,
00099 (int) (Y->stride));
00100 return BPM_SUCCESS;
00101 }
00102 else
00103 {
00104 bpm_error("invalid length in gsl_blas_dgemv(..)",
00105 __FILE__, __LINE__);
00106 return BPM_FAILURE;
00107 }
00108 }
00109
00110 void
00111 cblas_dgemv (const enum CBLAS_ORDER order, const enum CBLAS_TRANSPOSE TransA,
00112 const int M, const int N, const double alpha, const double *A,
00113 const int lda, const double *X, const int incX,
00114 const double beta, double *Y, const int incY)
00115 {
00116 int i, j;
00117 int lenX, lenY;
00118
00119 const int Trans = (TransA != CblasConjTrans) ? TransA : CblasTrans;
00120
00121 if (M == 0 || N == 0)
00122 return;
00123
00124 if (alpha == 0.0 && beta == 1.0)
00125 return;
00126
00127 if (Trans == CblasNoTrans) {
00128 lenX = N;
00129 lenY = M;
00130 } else {
00131 lenX = M;
00132 lenY = N;
00133 }
00134
00135
00136 if (beta == 0.0) {
00137 int iy = OFFSET(lenY, incY);
00138 for (i = 0; i < lenY; i++) {
00139 Y[iy] = 0.0;
00140 iy += incY;
00141 }
00142 } else if (beta != 1.0) {
00143 int iy = OFFSET(lenY, incY);
00144 for (i = 0; i < lenY; i++) {
00145 Y[iy] *= beta;
00146 iy += incY;
00147 }
00148 }
00149
00150 if (alpha == 0.0)
00151 return;
00152
00153 if ((order == CblasRowMajor && Trans == CblasNoTrans)
00154 || (order == CblasColMajor && Trans == CblasTrans)) {
00155
00156 int iy = OFFSET(lenY, incY);
00157 for (i = 0; i < lenY; i++) {
00158 double temp = 0.0;
00159 int ix = OFFSET(lenX, incX);
00160 for (j = 0; j < lenX; j++) {
00161 temp += X[ix] * A[lda * i + j];
00162 ix += incX;
00163 }
00164 Y[iy] += alpha * temp;
00165 iy += incY;
00166 }
00167 } else if ((order == CblasRowMajor && Trans == CblasTrans)
00168 || (order == CblasColMajor && Trans == CblasNoTrans)) {
00169
00170 int ix = OFFSET(lenX, incX);
00171 for (j = 0; j < lenX; j++) {
00172 const double temp = alpha * X[ix];
00173 if (temp != 0.0) {
00174 int iy = OFFSET(lenY, incY);
00175 for (i = 0; i < lenY; i++) {
00176 Y[iy] += temp * A[lda * j + i];
00177 iy += incY;
00178 }
00179 }
00180 ix += incX;
00181 }
00182 } else {
00183 bpm_error("unrecognised operation in cblas_dgemv(..)",
00184 __FILE__, __LINE__);
00185 return;
00186 }
00187 }