00001 #ifndef MATRIX_MUTILS_H
00002 #define MATRIX_MUTILS_H
00003
00004 #ifdef __cplusplus
00005 extern "C" {
00006 #endif
00007
00008 #include <Rdefines.h>
00009 #include <Rconfig.h>
00010 #include <R.h>
00011
00012 #ifdef ENABLE_NLS
00013 #include <libintl.h>
00014 #define _(String) dgettext ("Matrix", String)
00015 #else
00016 #define _(String) (String)
00017 #endif
00018
00019
00020 enum CBLAS_ORDER {CblasRowMajor=101, CblasColMajor=102};
00021 enum CBLAS_TRANSPOSE {CblasNoTrans=111, CblasTrans=112, CblasConjTrans=113};
00022 enum CBLAS_UPLO {CblasUpper=121, CblasLower=122};
00023 enum CBLAS_DIAG {CblasNonUnit=131, CblasUnit=132};
00024 enum CBLAS_SIDE {CblasLeft=141, CblasRight=142};
00025 #define RMJ CblasRowMajor
00026 #define CMJ CblasColMajor
00027 #define NTR CblasNoTrans
00028 #define TRN CblasTrans
00029 #define CTR CblasConjTrans
00030 #define UPP CblasUpper
00031 #define LOW CblasLower
00032 #define NUN CblasNonUnit
00033 #define UNT CblasUnit
00034 #define LFT CblasLeft
00035 #define RGT CblasRight
00036
00037 char norm_type(char *typstr);
00038 char rcond_type(char *typstr);
00039 double get_double_by_name(SEXP obj, char *nm);
00040 SEXP set_double_by_name(SEXP obj, double val, char *nm);
00041 SEXP as_det_obj(double val, int log, int sign);
00042 SEXP get_factors(SEXP obj, char *nm);
00043 SEXP set_factors(SEXP obj, SEXP val, char *nm);
00044 SEXP dgCMatrix_set_Dim(SEXP x, int nrow);
00045 int csc_unsorted_columns(int ncol, const int p[], const int i[]);
00046 void csc_sort_columns(int ncol, const int p[], int i[], double x[]);
00047 SEXP triple_as_SEXP(int nrow, int ncol, int nz,
00048 const int Ti [], const int Tj [], const double Tx [],
00049 char *Rclass);
00050 SEXP csc_check_column_sorting(SEXP A);
00051 void csc_compTr(int m, int n, int nnz,
00052 const int xp[], const int xi[], const double xx[],
00053 int ap[], int ai[], double ax[]);
00054 void ssc_symbolic_permute(int n, int upper, const int perm[],
00055 int Ap[], int Ai[]);
00056 double *nlme_symmetrize(double *a, const int nc);
00057 void nlme_check_Lapack_error(int info, const char *laName);
00058 SEXP nlme_replaceSlot(SEXP obj, SEXP names, SEXP value);
00059 SEXP nlme_weight_matrix_list(SEXP MLin, SEXP wts, SEXP adjst, SEXP MLout);
00060 SEXP Matrix_make_named(int TYP, char **names);
00061 SEXP check_scalar_string(SEXP sP, char *vals, char *nm);
00062 double *packed_to_full(double *dest, const double *src, int n,
00063 enum CBLAS_UPLO uplo);
00064 double *full_to_packed(double *dest, const double *src, int n,
00065 enum CBLAS_UPLO uplo, enum CBLAS_DIAG diag);
00066 double *packed_getDiag(double *dest, SEXP x);
00067
00068 extern
00069 #include "Syms.h"
00070
00071
00072 #define AZERO(x, n) {int _I_, _SZ_ = (n); for(_I_ = 0; _I_ < _SZ_; _I_++) (x)[_I_] = 0;}
00073
00074
00075 #define PACKED_LENGTH(n) ((n) * ((n) + 1))/2
00076
00085 static R_INLINE
00086 int packed_ncol(int len)
00087 {
00088 int disc = 8 * len + 1;
00089 int sqrtd = (int) sqrt((double) disc);
00090
00091 if (len < 0 || disc != sqrtd * sqrtd)
00092 error(_("invalid 'len' = %d in packed_ncol"));
00093 return (sqrtd - 1)/2;
00094 }
00095
00111 static R_INLINE
00112 SEXP ALLOC_SLOT(SEXP obj, SEXP nm, SEXPTYPE type, int length)
00113 {
00114 SEXP val = allocVector(type, length);
00115
00116 SET_SLOT(obj, nm, val);
00117 return val;
00118 }
00119
00130 static R_INLINE
00131 int* expand_cmprPt(int ncol, const int mp[], int mj[])
00132 {
00133 int j;
00134 for (j = 0; j < ncol; j++) {
00135 int j2 = mp[j+1], jj;
00136 for (jj = mp[j]; jj < j2; jj++) mj[jj] = j;
00137 }
00138 return mj;
00139 }
00140
00141
00157 static R_INLINE int
00158 check_csc_index(const int p[], const int i[], int row, int col, int missing)
00159 {
00160 int k, k2 = p[col + 1];
00161
00162 for (k = p[col]; k < k2; k++) if (i[k] == row) return k;
00163 if (!missing)
00164 error("row %d and column %d not defined in rowind and colptr",
00165 row, col);
00166 return missing;
00167 }
00168
00169 SEXP alloc3Darray(SEXPTYPE mode, int nrow, int ncol, int nface);
00170
00180 static R_INLINE
00181 int Lind(int k, int i)
00182 {
00183 if (k < i) error("Lind(k = %d, i = %d) must have k >= i", k, i);
00184 return (k * (k + 1))/2 + i;
00185 }
00186
00195 static R_INLINE
00196 int match_mat_dims(const int xd[], const int yd[])
00197 {
00198 return xd[0] == yd[0] && xd[1] == yd[1];
00199 }
00200
00201 double *expand_csc_column(double *dest, int m, int j,
00202 const int Ap[], const int Ai[], const double Ax[]);
00203
00211 static R_INLINE void
00212 int_permute(int i[], int n, const int perm[])
00213 {
00214 int j;
00215 for (j = 0; j < n; j++) i[j] = perm[i[j]];
00216 }
00217
00225 static R_INLINE void
00226 make_upper_triangular(int i[], int j[], int nnz)
00227 {
00228 int k;
00229 for (k = 0; k < nnz; k++) {
00230 if (i[k] > j[k]) {
00231 int tmp = i[k];
00232 i[k] = j[k];
00233 j[k] = tmp;
00234 }
00235 }
00236 }
00237
00238 void make_array_triangular(double *x, SEXP from);
00239
00240 #ifdef __cplusplus
00241 }
00242 #endif
00243
00244 #endif /* MATRIX_MUTILS_H_ */