Арифметика высокой точности в ALGLIB, как ее использовать? Опции

[MSP430F]

Всем доброго времени суток. Надеюсь, что мой вопрос не останется без ответа.
Использую функцию rmatrixsolvels из пакета alglib. Подключил просто - закинул в папку с проектом все хедеры и сишники alglib. Из 13 файлов пришлось подключить 6 для вызова лишь одной функции. exe подрос заметно в размерах, время компиляции увеличилось заметно, ну да ладно, главное цель достигнута. Решается система линейных уравнений с прямоугольной матрицей - число уравнений превышает число неизвестных. Минимизируется невязка. Но чувствую, не хватает точности. В alglib есть возможность работы с арифметикой высокой точности http://alglib.sources.ru/equations/linear.php Но как ее подключить ? Что-то я не соображу... Так все запутано в этом alglib. Использую компилятор MinGW-W64. Буду признателен за любые подсказки.

[thermit]

MSP430F:
Но чувствую, не хватает точности.

Каким образом почувствовали? Какого размера система? Каким способом решается?
Вероятно, необходимость в высокой точности необоснована.

[Guest_TSerg_*]

Проверьте с использованием мат.пакетов Матлаб, Маткад..

[MSP430F]

Каким образом почувствовали? Какого размера система? Каким способом решается?
Вероятно, необходимость в высокой точности необоснована.

Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

[Xenia]

Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

Если у вас на квадратной матрице rmatrixsolvels дает ошибку от 10-5 до 10-3, то при переходе на прямоугольную матрицу хуже уже не будет. Ухудшение на "3-5 порядков" - это плата за метод более общего типа (в нем арифметики больше), которую вы уже заплатили.

[thermit]

MSP430F:
Я проверил максимальную ошибку в вычислении свободных членов по найденным решениям, использую функции rmatrixsolve (для квадратной матрицы) и rmatrixsolvels (также в режиме квадратной матрицы). Ошибки отличаются на 3-5 порядков (в первой от 10-10 до 10-7, во второй от 10-5 до 10-3), хотя данные идентичны. Мне же предстоит решать систему не с квадратной, а прямоугольной матрицей. Если с квадратной такие ошибки,
то что же ожидать от прямоугольной ?

Исходные данные можете предоставить?

[MSP430F]

Если у вас на квадратной матрице rmatrixsolvels дает ошибку от 10-5 до 10-3, то при переходе на прямоугольною матрицу хуже уже не будет. Ухудшение на "3-5 порядков" - это плата за метод более общего типа (в нем арифметики больше), которую вы уже заплатили.

Спасибо за ответ! Я чувствовал нечто подобное, а Вы подтвердили мои догадки. Да, похоже, что использовать арифметику высокой точности по-любому очень проблемно здесь. Но если бы кто, подсказал, как, то можно было и попробовать. Если бы можно было просто заменить double на long double... Но там иакой финт не пройдет, слишком все в alglib наворочено.

Исходные данные можете предоставить?

Да. В файле матрица 21x21 элемент и массив B из 21-го элемента (правая часть). Максимальная ошибка для этих данных соответственно у меня получилась 1,1х10-9 и 1,9х10-4.

Прикрепленные файлы

 Matrix_A_and_array_B.txt ( 4.13 килобайт ) Кол-во скачиваний: 44

 

[thermit]

1. Плохо обусловленная задача. Число обусловленности матрицы ~10^8. При погрешности входных данных 10^-16 ошибка в результате <=~10^-8
2. Боюсь, алглиб вас не спасет. Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. На вашем примере гивенс дает максимальную погрешность 1.5*10^-9.

ps
Ну или ковыряться с расширенной точностью, если ошибка 10^-8 вас не устраивает.

[Xenia]

Боюсь, алглиб вас не спасет. Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой. На вашем примере гивенс дает максимальную погрешность 1.5*10^-9.

Alglib действительно не спасет, там алгоритм для этого дела гадкий . За Хаусхолдера же бояться нечего, это ортогональное преобразование, и ошибок оно практически не вносит (лишь на уровне точности арифметики). А вот на стадии обратной подстановки при плохо обусловленной матрице неминуемо возникнут проблемы, т.к. обусловленность полученной треугольной матрицы останется такой же малой, как у исходной.

Боюсь, что в данном случае надо считать через сингулярное разложение исходной матрицы, а то и понижая при этом размерность. Т.е. получать псевдорешение.

Впрочем, я бы рекомендовала сначала не заморачиваться всем этим, а применить процедуру rmatrixsolvels, как она есть, к прямоугольной матрице, на которой, по моим ощущениям, ошибка должна ученьшиться по сравнению с квадратной.

[MSP430F]

Вам нужно свести исходную матрицу к треугольной при помощи последовательности отражений (преобразование хаусхолдера) или плоских вращений (преобразование гивенса) и найти решение обратной подстановкой.

И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной...

Боюсь, что в данном случае надо считать через сингулярное разложение исходной матрицы.

Но ведь Alglib так и считает...

Сообщение отредактировал MSP430F - Dec 10 2013, 16:45

[thermit]

MSP430F:
И где же взять код этих адгоритмов для прямоугольной матрицы, не подскажите? Для квадратной алгоритм вращений у меня есть, он очень неплох. А вот для прямоугольной...

Можно легко модифицировать для прямоугольной. Если исходники есть.

ps

Могу нарисовать за умеренную плату, гы...

[MSP430F]

Можно легко модифицировать для прямоугольной. Если исходники есть.

Вот гляньте этот файл. Эта функция дает ошибку 2.7*10-8 на тех же исходных данных. Она не оригинальна, чуть модифицирована отсюда Мной добавлено создание копии матрицы в теле функции, а создание массива результатов убрано. Также введена проверка деления на 0.

Сообщение отредактировал MSP430F - Dec 11 2013, 07:54

Прикрепленные файлы

 Rotation.txt ( 1.74 килобайт ) Кол-во скачиваний: 40

 

[thermit]

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:
Ошибка ~9*10^-10

main.c

Код

#include <stdio.h>
#include <math.h>

int SolveByGivensRotation(double *A,double *b,int R,int C,double *x);





double a[21][21]=
{
{1.000000,    1.229402,    1.229566,    1.230016,    1.229648,    1.229771,    1.229975,    1.229361,    1.229525,    1.229648,    1.229484,    1.229238,    1.229156,    1.228992,    1.228214,    1.227518,    1.227928,    1.227108,    1.227231,    1.226453,    1.226576},    
{1.000000,    1.229220,    1.228551,    1.229852,    1.229774,    1.230227,    1.230539,    1.230034,    1.229774,    1.229279,    1.228967,    1.229019,    1.228978,    1.228677,    1.228413,    1.227573,    1.227919,    1.227521,    1.226982,    1.226752,    1.226759},    
{1.000000,    1.229566,    1.229689,    1.229846,    1.230118,    1.229974,    1.229979,    1.229705,    1.229582,    1.229385,    1.229313,    1.229417,    1.228537,    1.228060,    1.228113,    1.227828,    1.227652,    1.227244,    1.227276,    1.226674,    1.226855},    
{1.000000,    1.228701,    1.229879,    1.229808,    1.230219,    1.230501,    1.230044,    1.229808,    1.229361,    1.229078,    1.229125,    1.229088,    1.228816,    1.228577,    1.227816,    1.228129,    1.227769,    1.227281,    1.227072,    1.227079,    1.226659},    
{1.000000,    1.229696,    1.229674,    1.229618,    1.229730,    1.229786,    1.229718,    1.229696,    1.229607,    1.229741,    1.229528,    1.229237,    1.229349,    1.228286,    1.227962,    1.227481,    1.227738,    1.227279,    1.226708,    1.226932,    1.227044},    
{1.000000,    1.229095,    1.229253,    1.230450,    1.229750,    1.229885,    1.229389,    1.229389,    1.229479,    1.229614,    1.229547,    1.229140,    1.229253,    1.228824,    1.228576,    1.227943,    1.227966,    1.227153,    1.226860,    1.226815,    1.226431},    
{1.000000,    1.228654,    1.228773,    1.230379,    1.230498,    1.229963,    1.229784,    1.229190,    1.230022,    1.229784,    1.229071,    1.228773,    1.228952,    1.228297,    1.228714,    1.228119,    1.227822,    1.226810,    1.227286,    1.226929,    1.226989},    
{1.000000,    1.229305,    1.229601,    1.229552,    1.229848,    1.230292,    1.229453,    1.229552,    1.229700,    1.229798,    1.229650,    1.228713,    1.228368,    1.229108,    1.228516,    1.227480,    1.228072,    1.227185,    1.227037,    1.227234,    1.226346},    
{1.000000,    1.228883,    1.228629,    1.230443,    1.230170,    1.230034,    1.230053,    1.230190,    1.229663,    1.228766,    1.228766,    1.228961,    1.229000,    1.228746,    1.228473,    1.227440,    1.227849,    1.228122,    1.227089,    1.226640,    1.226893},    
{1.000000,    1.229543,    1.229630,    1.229356,    1.229499,    1.229683,    1.230001,    1.229833,    1.229839,    1.229518,    1.229374,    1.229143,    1.229059,    1.229181,    1.228151,    1.227593,    1.227655,    1.227322,    1.227116,    1.226638,    1.226676},    
{1.000000,    1.229429,    1.229542,    1.230312,    1.229936,    1.230049,    1.230105,    1.229336,    1.229148,    1.229054,    1.228716,    1.228491,    1.229636,    1.228359,    1.228453,    1.227590,    1.228266,    1.228040,    1.227609,    1.226276,    1.226464},    
{1.000000,    1.229692,    1.229447,    1.229575,    1.229739,    1.230022,    1.229872,    1.229878,    1.229592,    1.229464,    1.229258,    1.229183,    1.229292,    1.228375,    1.227878,    1.227933,    1.227636,    1.227453,    1.227027,    1.227061,    1.226433},    
{1.000000,    1.229186,    1.229563,    1.229939,    1.230316,    1.229971,    1.229469,    1.229312,    1.229374,    1.229531,    1.229469,    1.229280,    1.229061,    1.228527,    1.228810,    1.227617,    1.227335,    1.226644,    1.226833,    1.227021,    1.227554},    
{1.000000,    1.229820,    1.229803,    1.230004,    1.229990,    1.229898,    1.229716,    1.229279,    1.229198,    1.229403,    1.229509,    1.229152,    1.229086,    1.228881,    1.227502,    1.227404,    1.227758,    1.227441,    1.227067,    1.226886,    1.227012},    
{1.000000,    1.229545,    1.229452,    1.229675,    1.229582,    1.229861,    1.229805,    1.229712,    1.229694,    1.229824,    1.229285,    1.229174,    1.229025,    1.228413,    1.228394,    1.228227,    1.227707,    1.227206,    1.227002,    1.226575,    1.226649},    
{1.000000,    1.229498,    1.229330,    1.229305,    1.229378,    1.229955,    1.230075,    1.229882,    1.229642,    1.229522,    1.229618,    1.229354,    1.228993,    1.228176,    1.228945,    1.227743,    1.227310,    1.227671,    1.226949,    1.226829,    1.226637},    
{1.000000,    1.229748,    1.229748,    1.229999,    1.229971,    1.229776,    1.229608,    1.229385,    1.229469,    1.229497,    1.229329,    1.229134,    1.228799,    1.228687,    1.228576,    1.227516,    1.227237,    1.227432,    1.227460,    1.227153,    1.226288},    
{1.000000,    1.228414,    1.229260,    1.229719,    1.229719,    1.230036,    1.230000,    1.229754,    1.229542,    1.229260,    1.229366,    1.229401,    1.229190,    1.228943,    1.228520,    1.228379,    1.228238,    1.226899,    1.226547,    1.226793,    1.226828},    
{1.000000,    1.229178,    1.228977,    1.229153,    1.230486,    1.229706,    1.229857,    1.229304,    1.229304,    1.229405,    1.229556,    1.229480,    1.229027,    1.229153,    1.228675,    1.228398,    1.227694,    1.227719,    1.226814,    1.226487,    1.226436},    
{1.000000,    1.229249,    1.229674,    1.229702,    1.229759,    1.229957,    1.229759,    1.229702,    1.229787,    1.229702,    1.229475,    1.229192,    1.228399,    1.227946,    1.227776,    1.227890,    1.228314,    1.227691,    1.228059,    1.226360,    1.226417},    
{1.000000,    1.228960,    1.229996,    1.230061,    1.230319,    1.230190,    1.230190,    1.229219,    1.229414,    1.228960,    1.228766,    1.228572,    1.229090,    1.229155,    1.227990,    1.227860,    1.227666,    1.227278,    1.226890,    1.227278,    1.226955}    
};

double b[21]=
{
1.226985,    1.226294,    1.226775,    1.225659,    1.227156,    1.226476,    1.225978,    1.225508,    1.226484,    1.225971,    1.226426,    1.226622,    1.224824,    1.226897,    1.226687,    1.226372,    1.226260,    1.226441,    1.226009,    1.226643,    1.226760
};

double aorig[21][21],aa[21][21],x[21],bb[21];

int main()
{
    int i,j;
    double acc,maxerror;

    for(i=0;i<21;i++)
    {
        for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];
        bb[i]=b[i];
    }


SolveByGivensRotation(&aa[0][0],bb,21,21,x);


maxerror=0.0;
for(i=0;i<21;i++)
{
     acc=b[i];
     for(j=0;j<21;j++) acc-=aorig[i][j]*x[j];
     if(fabs(acc)>maxerror) maxerror=fabs(acc);
}

printf("Max error  = %.15g\n",maxerror);

return 0;

}

givens.c

Код

#include <math.h>

#define SQR(x)     (x*x)




int Rotation(double *A,double *b,int R,int C,int ri,int ci)
{
    int i;
    double t1,t2,t3,c,s;

    t1=A[C*ci+ci];
    t2=A[C*ri+ci];

    t3=sqrt(SQR(t1)+SQR(t2));
    if(t3==0.0) return 1;
    c=t1/t3;
    s=t2/t3;
    for(i=0;i<C;i++)
    {
        t1=c*A[ci*C+i]+s*A[ri*C+i];
        t2=c*A[ri*C+i]-s*A[ci*C+i];
        A[ci*C+i]=t1;
        A[ri*C+i]=t2;
    }
    t1=c*b[ci]+s*b[ri];
    t2=c*b[ri]-s*b[ci];
    b[ci]=t1;
    b[ri]=t2;

    return 0;
}




int GivensTriang(double *A,double *b,int R,int C)
{
    int i,j;


    for(i=0;i<C;i++)//For col
        for(j=R-1;j>i;j--) //For row
            if(Rotation(A,b,R,C,j,i)) return 1;

    return 0;
}

void BackwardSubstitution(double *A,double *b,double *x,int R,int C)
{
int i,j;
double tmp;


if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1];
for(i=C-2;i>=0;i--)
{
     tmp=b[i];
     for(j=i+1;j<C;j++) tmp-=x[j]*A[i*C+j];
     if(A[i*C+i]==0.0) x[i]=0.0;
     else x[i]=tmp/A[i*C+i];
}
}

int SolveByGivensRotation(double *A,double *b,int R,int C,double *x)
{
    if(GivensTriang(A,b,R,C)) return 1;
    BackwardSubstitution(A,b,x,R,C);

    return 0;
}

Сообщение отредактировал thermit - Dec 11 2013, 12:48

[_pv]

Код

void GaussSolve(double * A, double * Y, int num){
  //zeroeing left lower part of A[]
  for (int k = 0; k < num - 1; k++){   //diagonal
    for (int j = k + 1; j < num; j++){  //vertical
      double m = A[k + j * num] / A[k + k * num];
      for (int i = k+1; i < num; i++)  A[i + j * num] -= A[i + k * num] * m;  //horizontal
      Y[j] -= Y[k] * m;
    }
  }
  //zeroeing right upper part of A[], Y only.
  for (int k = num - 1; k > 0; k--)
    for (int j = k - 1; j >= 0; j--)
      Y[j] -= Y[k] * A[k + j * num] / A[k + k * num];
  //return answer in Y
  for (int i = 0; i < num; i++) Y[i] /= A[i + i * num];
}

обычный Гаусс на этих данных даёт ошибку 10^-12.
для прямоугольной матрицы при нахождении (A* A)^-1 A* поди не сильно больше наберётся.

[thermit]

_pv:
обычный Гаусс на этих данных даёт ошибку 10^-12.

Не дает. Те же самые ~9*10^-10.
Матрицу нужно транспонировать.

[MSP430F]

Честно говоря, лень ковыряться с чужими исходниками.
Должно быть как-то так:

Это Вы мне что, прямо таки код рабочий дали ? Спасибо большое, буду проверять. Если все OK, то куда donation перевести ?

Так как же должно быть в BackwardSubstitution():

if(A[(C-1)*R+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*R+C-1]; или

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1]; ?


Сообщение отредактировал MSP430F - Dec 11 2013, 13:27

[_pv]

5*10^-10
если то же самое в чём радость по сравнению с Гауссом?

[thermit]

Это Вы мне что, прямо таки код рабочий дали ? Спасибо большое, буду проверять. Если все OK, то куда donation перевести ?

Так как же должно быть в BackwardSubstitution():

if(A[(C-1)*R+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*R+C-1]; или

if(A[(C-1)*C+C-1]==0.0) x[C-1]=0.0;
else x[C-1]=b[C-1]/A[(C-1)*C+C-1]; ?

Последняя редакция верна. Первый раз скопировал с ошибкой.

_pv:
если то же самое в чём радость по сравнению с Гауссом?

На плохо обусловленных системах с числом уравнений большим числа неизвестных ошибка гаусса будет на порядок больше. Нужно будет принимать специальные меры.
Кроме того, применение ортогонального преобразования гарантирует состоятельность и численную устойчивость алгоритма. Про гаусса такого сказать нельзя.
Если число операций не является приоритетом, лучше использовать вращения или отражения.

зы
старик уилкинсон утверждает, что ему всегда хватало гаусса. Но это было в эпоху медленных компьютеров...

[Guest_TSerg_*]

старик уилкинсон утверждает, что ему всегда хватало гаусса. Но это было в эпоху медленных компьютеров...

Не только.
Все чаще дело имеем с разреженными системами, а Гаусс разреженность лучше сохраняет, в отличии от методов отражений и вращений.

[MSP430F]

Я правильно догадался, что в этом месте ошибка:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i];

Должно быть:

for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[i][j];

Первый индекс - строка, второй - столбец. Это же одинаково для всех массивов ?

[thermit]

Ошибки здесь нет. Исходная матрица лежит в памяти постолбцово (1-я строка состоит из единиц). Во всяком случае я так считал. Для работы алгоритмов матрица должна лежать в памяти построчно. for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i]; Это и есть такое преобразование (транспонирование).

Сообщение отредактировал thermit - Dec 12 2013, 09:19

[MSP430F]

Ошибки здесь нет. Исходная матрица лежит в памяти постолбцово (1-я строка состоит из единиц). Во всяком случае я так считал. Для работы алгоритмов матрица должна лежать в памяти построчно. for(j=0;j<21;j++) aa[i][j]=aorig[i][j]=a[j][i]; Это и есть такое преобразование (транспонирование).

Если считать, что первый индекс - строка, а второй - столбец, то первая строка как раз и не состоит из 1.0, а равна
{1.000000, 1.229402, 1.229566, 1.230016, 1.229648, 1.229771, 1.229975, 1.229361, 1.229525, 1.229648, 1.229484, 1.229238, 1.229156, 1.228992, 1.228214, 1.227518, 1.227928, 1.227108, 1.227231, 1.226453, 1.226576},
как она и проинициализирована. Зачем же здесь тогда транспонирование ?

Можете проверить:

i = 0;
for(j=0; j<21; j++) printf("a = %f\n", a[i][j]);

Сообщение отредактировал MSP430F - Dec 12 2013, 10:21

[_pv]

тогда простой Гаусс даст ошибку в 10^-12.
а для прямоугольной матрицы тем же Гауссом постойте псевдообратную и будет вам счастье и без ALGLIBа.

[thermit]

MSP430F:
Зачем же здесь тогда транспонирование ?

Для данных, которые вы прислали погрешность вычислялась (и с вашей совпадала) для матрицы, где строки являются столбцами из вашей таблицы. Поскольку я вставил ее себе в код копипастом, пришлось ввести транспонирование. С индексами все путем: 1-й - номер строки, 2-й - номер столбца.
Честно говоря, не понял проблемы. Разместите свои данные построчно да и все дела.

_pv:
будет вам счастье и без ALGLIBа.

Собственно, то, что алглиб тут лишний, давно понятно.

[MSP430F]

Для данных, которые вы прислали погрешность вычислялась (и с вашей совпадала) для матрицы, где строки являются столбцами из вашей таблицы. Поскольку я вставил ее себе в код копипастом, пришлось ввести транспонирование. С индексами все путем: 1-й - номер строки, 2-й - номер столбца.
Честно говоря, не понял проблемы. Разместите свои данные построчно да и все дела.

Нет, Вы не правы. В моих данных все на месте - строки это строки, столбцы - столбцы. Первый элемент каждой строки - это 1.0 (так и должно быть). Транспонирование не нужно.

[thermit]

MSP430F:
Нет, Вы не правы. В моих данных все на месте - строки это строки, столбцы - столбцы. Первый элемент каждой строки - это 1.0 (так и должно быть). Транспонирование не нужно.

Это непринципиально. Кстати, что за данные? Откуда?

[MSP430F]

Это непринципиально. Кстати, что за данные? Откуда?

Проверяю метод ближайших соседей для прогнозирования значения временного ряда хаотического типа. Для примера взял котировки EUR/USD с шагом 1 час.

Ваш код хорошо работает. Правда, пришлось немного попотеть, чтобы можно было в качестве аргумента передавать матрицу в виде нормального двумерного массива, так как размерность матрицы может меняться динамически. Ваш код дает ошибку не хуже, чем alglib. При небольшой прямоугольности разница заметна в пользу Вашего кода. Если матрица 21 * 25 или 21 * 31, то разница в ошибке незаметна. Зато код получается небольшой и нет кучи warning-ов при компиляции, как при использовании alglib.

[thermit]

MSP430F:
Правда, пришлось немного попотеть, чтобы можно было в качестве аргумента передавать матрицу в виде нормального двумерного массива

Зачем передавать в виде двумерного массива? Решателю достаточно передать адрес начала матрицы и ее размерность.

[Xenia]

Зачем передавать в виде двумерного массива? Решателю достаточно передать адрес начала матрицы и ее размерность.

У топик-стартера матрица аллокируется не одним блоком памяти, а построчно:

Код

  double **copy_mass = new double *[cnt_str]; // создаем копию матрицы, так как элементы матрицы будут модифицированы
  for(i = 0; i < cnt_str; i++){
    copy_mass[i] = new double [cnt_str + 1]; // N столбцов = N строк + 1
    for(j = 0; j < cnt_str + 1; j++) copy_mass[i][j] = mass[i][j];
  }

Однако это обстоятельство ничуть не мешает ее передавать через указатель double**, а еще лучше, как &double**.
Хотя на C++ можно было бы создать для матрицы класс, где кроме указателя double**, хранились бы и обе ее размерности. А конструктор и деструктор память под эту матрицу аллокировали и деаллокировали. В этом случае ссылку (здесь она удобнее указателя) на матрицу можно передавать одним единственным парамером. Кроме того, данный класс можно снабдить фукциями членами класса дял матричного умножения, сложения и т.п.