Вызов подпрограммы LAPACK непосредственно в RcppArmadillo - PullRequest
Новая
       43

Вызов подпрограммы LAPACK непосредственно в RcppArmadillo

0 голосов
/ 20 октября 2018

Поскольку у Armadillo (afaik) нет треугольного решателя, я бы хотел использовать треугольный решатель LAPACK, доступный в dtrtrs.Я посмотрел на следующие два ( first , second ) SO потока и сложил что-то вместе, но это не работает.

Я создал свежий пакетиспользуя RStudio, одновременно включив RcppArmadillo.У меня есть файл заголовка header.h: 

#include <RcppArmadillo.h>

#ifdef ARMA_USE_LAPACK
#if !defined(ARMA_BLAS_CAPITALS)
#define arma_dtrtrs dtrtrs
#else
#define arma_dtrtrs DTRTRS
#endif
#endif

extern "C" {
  void arma_fortran(arma_dtrtrs)(char* UPLO, char* TRANS, char* DIAG, int* N, int* NRHS,
                    double* A, int* LDA, double* B, int* LDB, int* INFO);
}

int trtrs(char uplo, char trans, char diag, int n, int nrhs, double* A, int lda, double* B, int ldb);

static int trisolve(const arma::mat &in_A, const arma::mat &in_b, arma::mat &out_x);

, который, по сути, является ответом на первый связанный вопрос, а также с функцией-оберткой и основной функцией.Основная часть функций указана в trisolve.cpp и выглядит следующим образом: 

#include "header.h"

int trtrs(char uplo, char trans, char diag, int n, int nrhs, double* A, int lda, double* B, int ldb) {
  int info = 0;
  wrapper_dtrtrs_(&uplo, &trans, &diag, &n, &nrhs, A, &lda, B, &ldb, &info);
  return info;
}


static int trisolve(const arma::mat &in_A, const arma::mat &in_b, arma::mat &out_x) {
  size_t  rows = in_A.n_rows;
  size_t  cols = in_A.n_cols;

  double *A = new double[rows*cols];
  double *b = new double[in_b.size()];

  //Lapack has column-major order
  for(size_t col=0, D1_idx=0; col<cols; ++col)
  {
    for(size_t row = 0; row<rows; ++row)
    {
      // Lapack uses column major format
      A[D1_idx++] = in_A(row, col);
    }
    b[col] = in_b(col);
  }

  for(size_t row = 0; row<rows; ++row)
  {
    b[row] = in_b(row);
  }

  int info = trtrs('U', 'N', 'N', cols, 1, A, rows, b, rows);

  for(size_t col=0; col<cols; col++) {
    out_x(col)=b[col];
  }

  delete[] A;
  delete[] b;

  return 0;
}


// [[Rcpp::export]]

arma::mat RtoRcpp(arma::mat A, arma::mat b) {
  arma::uword n = A.n_rows;
  arma::mat x = arma::mat(n, 1, arma::fill::zeros);

  int info = trisolve(A, b, x);
  return x;
}

У меня есть (по крайней мере) две проблемы:

  1. При попытке компиляции яполучить: conflicting types for 'dtrtrs_' из заголовочного файла.Однако я не вижу, что не так с входными данными (и это буквально скопировано из второго связанного потока).
  2. Неудивительно, что wrapper_dtrtrts_ не правильно.Но из того, что я могу сказать из compiler_setup.hpp, Армадилло должен создать для меня функцию с именем wrapper_dtrtrs_.Какое имя я должен использовать здесь в основном cpp файле?

Ответы [ 2 ]

0 голосов
/ 21 октября 2018

Armadillo уже использует dtrtrs для решения трехдиагональных задач.Некоторые ссылки на код:

Так что, если мы сможем запустить этот оператор отладки, мы можем быть уверены, что dtrtrs действительно используется: 

#define ARMA_EXTRA_DEBUG
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

// [[Rcpp::export]]
void testTrisolve() {
  arma::mat A = arma::randu<arma::mat>(5,5);
  arma::mat B = arma::randu<arma::mat>(5,5);

  arma::mat X1 = arma::solve(A, B);
  arma::mat X3 = arma::solve(arma::trimatu(A), B);
}

/*** R
testTrisolve()
*/

Это создает много отладочных сообщений, среди них:

lapack::gesvx()
[...]
lapack::trtrs()

Итак, мы ясно видим, что в трехдиагональном случае используется dtrtrs.

Что касается ваших исходных вопросов:

  1. Ошибка конфликтующих типов:следствие того, что Арамдильо уже использует dtrtrs, но с немного другой сигнатурой (A - const).
  2. Имя уровня C для функции Fortran зависит от значений ARMA_BLAS_UNDERSCORE и ARMA_USE_WRAPPER.Я не уверен, так ли это всегда, но для меня первое определено, а второе нет (ср. config.hpp), что приводит к dtrtrs_ в качестве имени.

Действительно, если ядобавьте const, где Armadillo использует его, и вызовите функцию как dtrtrs_, ваш код компилируется без ошибок или предупреждений (за исключением неиспользуемой переменной ...): 

// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>

extern "C" {
  void arma_fortran(dtrtrs)(char* UPLO, char* TRANS, char* DIAG, int* N, int* NRHS,
                    const double* A, int* LDA, double* B, int* LDB, int* INFO);
}

int trtrs(char uplo, char trans, char diag, int n, int nrhs, double* A, int lda, double* B, int ldb) {
  int info = 0;
  dtrtrs_(&uplo, &trans, &diag, &n, &nrhs, A, &lda, B, &ldb, &info);
  return info;
}

[...]
0 голосов
/ 21 октября 2018

У броненосца уже есть треугольный решатель.Код, адаптированный из документации : 

mat A(5,5, fill::randu);
// ... make A triangular here ...

mat B(5,5, fill::randu);

// tell solve() to treat A as upper triangular
// and automatically enable fast mode
mat X = solve(trimatu(A), B);

Исходя из документации, похоже, что решатель Armadillo может автоматически обнаруживать матрицы полос и симметричные положительно определенные матрицы.

...