Décomposition de LU à partir de recettes numériques ne fonctionne pas; Qu’est-ce que je fais mal?

J’ai littéralement copié et collé à partir du code source fourni pour Numerical Recipes for C pour la décomposition de masortingce LU sur place, le problème est que cela ne fonctionne pas.

Je suis sûr que je fais quelque chose de stupide mais apprécierais que quiconque puisse me diriger dans la bonne direction à ce sujet; J’ai travaillé toute la journée et je ne vois pas ce que je fais mal.

MISE À JOUR POST-ANSWER: Le projet est terminé et fonctionne . Merci à tous pour leurs conseils.

#include  #include  #include  #define MAT1 3 #define TINY 1e-20 int h_NR_LU_decomp(float *a, int *indx){ //Taken from Numerical Recipies for C int i,imax,j,k; float big,dum,sum,temp; int n=MAT1; float vv[MAT1]; int d=1.0; //Loop over rows to get implicit scaling info for (i=0;i<n;i++) { big=0.0; for (j=0;j big) big=temp; if (big == 0.0) return -1; //Singular Masortingx vv[i]=1.0/big; } //Outer kij loop for (j=0;j<n;j++) { for (i=0;i<j;i++) { sum=a[i*MAT1+j]; for (k=0;k<i;k++) sum -= a[i*MAT1+k]*a[k*MAT1+j]; a[i*MAT1+j]=sum; } big=0.0; //search for largest pivot for (i=j;i<n;i++) { sum=a[i*MAT1+j]; for (k=0;k= big) { big=dum; imax=i; } } //Do we need to swap any rows? if (j != imax) { for (k=0;k<n;k++) { dum=a[imax*MAT1+k]; a[imax*MAT1+k]=a[j*MAT1+k]; a[j*MAT1+k]=dum; } d = -d; vv[imax]=vv[j]; } indx[j]=imax; if (a[j*MAT1+j] == 0.0) a[j*MAT1+j]=TINY; for (k=j+1;k<n;k++) { dum=1.0/(a[j*MAT1+j]); for (i=j+1;i<n;i++) a[i*MAT1+j] *= dum; } } return 0; } void main(){ //3x3 Matrix float exampleA[]={1,3,-2,3,5,6,2,4,3}; //pivot array (not used currently) int* h_pivot = (int *)malloc(sizeof(int)*MAT1); int retval = h_NR_LU_decomp(&exampleA[0],h_pivot); for (unsigned int i=0; i<3; i++){ printf("\n%d:",h_pivot[i]); for (unsigned int j=0;j<3; j++){ printf("%.1lf,",exampleA[i*3+j]); } } } 

WolframAlpha dit que la réponse devrait être

 1,3,-2 2,-2,7 3,2,-2 

Je suis en train:

 2,4,3 0.2,2,-2.8 0.8,1,6.5 

Et jusqu’à présent, j’ai trouvé au moins 3 versions différentes du “même” algorithme, je suis donc complètement confus.

PS oui je sais qu’il y a au moins une douzaine de bibliothèques différentes pour le faire, mais je suis plus intéressé à comprendre ce que je fais mal que la bonne réponse.

PPS puisque, dans la décomposition de LU, la masortingce résultante inférieure est l’unité et qu’en utilisant l’algorithme de Crouts comme implémenté (je pense) implémenté, l’access aux index de tableau est toujours sécurisé, L et U peuvent être superposés l’un sur l’autre sur place; d’où la seule masortingce résultante pour cela.

    Je pense qu’il y a quelque chose de fondamentalement faux dans vos indices. Ils ont parfois des valeurs de début et de fin inhabituelles, et la boucle externe sur j au lieu de i me rend méfiant.

    Avant de demander à quelqu’un d’examiner votre code, voici quelques suggestions:

    • vérifiez vos index
    • se débarrasser de ces tentatives d’obfuscation en utilisant sum
    • utiliser une macro a(i,j) au lieu de a[i*MAT1+j]
    • écrire des sous-fonctions au lieu de commentaires
    • supprimer les parties inutiles, en isolant le code erroné

    Voici une version qui suit ces suggestions:

     #define MAT1 3 #define a(i,j) a[(i)*MAT1+(j)] int h_NR_LU_decomp(float *a, int *indx) { int i, j, k; int n = MAT1; for (i = 0; i < n; i++) { // compute R for (j = i; j < n; j++) for (k = 0; k < i-2; k++) a(i,j) -= a(i,k) * a(k,j); // compute L for (j = i+1; j < n; j++) for (k = 0; k < i-2; k++) a(j,i) -= a(j,k) * a(k,i); } return 0; } 

    Ses principaux avantages sont:

    • c'est lisible
    • Ça marche

    Cela manque de pivot, cependant. Ajoutez des sous-fonctions au besoin.

    Mon conseil: ne copiez pas le code de quelqu'un d'autre sans le comprendre.

    La plupart des programmeurs sont de mauvais programmeurs.

    Pour l’amour de tout ce qui est saint, n’utilisez pas de code Numerical Recipies sauf pour une implémentation de jouet à des fins pédagogiques des algorithmes décrits dans le texte – et, en réalité, le texte n’est pas si génial. Et, pendant que vous apprenez, le code non plus.

    Évidemment , ne mettez pas de routine Numerical Recipies dans votre propre code – la licence est incroyablement ressortingctive, en particulier en raison de la qualité du code. Vous ne pourrez pas dissortingbuer votre propre code si vous avez des éléments NR dedans.

    Vérifiez si une bibliothèque LAPACK est déjà installée sur votre système. C’est l’interface standard avec les routines d’algèbre linéaire en science et ingénierie computationnelles. Même si ce n’est pas parfait, vous pourrez trouver des bibliothèques lapack pour n’importe quelle machine sur laquelle vous déplacerez votre code, et vous pouvez simplement comstackr, lier et exécuter. S’il n’est pas déjà installé sur votre système, votre gestionnaire de paquets (rpm, apt-get, fink, port, etc.) connaît probablement lapack et peut l’installer à votre place. Sinon, tant que vous avez un compilateur Fortran sur votre système, vous pouvez le télécharger et le comstackr à partir d’ ici , et les liaisons en C standard peuvent être trouvées juste en dessous de la même page.

    La raison pour laquelle il est si pratique d’avoir une API standard pour les routines d’algèbre linéaire est qu’elles sont si courantes, mais leurs performances dépendent tellement du système. Ainsi, par exemple, Goto BLAS est une implémentation incroyablement rapide pour les systèmes x86 des opérations de bas niveau nécessaires à l’algèbre linéaire; Une fois que LAPACK fonctionne, vous pouvez installer cette bibliothèque pour que tout soit aussi rapide que possible.

    Une fois que vous avez installé LAPACK, la routine pour effectuer la factorisation LU d’une masortingce générale est SGETRF pour les floats ou DGETRF pour les doubles. Il existe d’autres routines plus rapides si vous connaissez la structure de la masortingce – qu’il s’agisse d’une définition positive symésortingque, par exemple (SBPTRF), ou d’une configuration sortingdiagonale (STDTRF). C’est une grande bibliothèque, mais une fois que vous avez appris à vous en servir, vous disposerez d’un équipement très puissant dans votre boîte à outils numériques.

    Ce qui me semble le plus méfiant est la partie intitulée “recherche du plus grand pivot”. Cela ne fait pas que chercher mais cela change aussi la masortingce A. J’ai du mal à croire que c’est correct.

    Les différentes versions de l’algorithme LU diffèrent en ce qui concerne le pivotement, assurez-vous donc de bien comprendre. Vous ne pouvez pas comparer les résultats de différents algorithmes. Une meilleure vérification consiste à vérifier si L fois U est égal à votre masortingce d’origine, ou une permutation de celle-ci si votre algorithme est pivotant. Cela étant dit, votre résultat est faux car le déterminant est faux (un pivot ne modifie pas le déterminant, à l’exception du signe).

    En dehors de cela, @Philip a de bons conseils. Si vous voulez comprendre le code, commencez par comprendre la décomposition de la LU sans pivoter.

    Pour paraphraser Albert Einstein:

    … Un homme avec une montre sait toujours l’heure exacte, mais un homme avec deux n’est jamais sûr ….

    Votre code ne produit définitivement pas le bon résultat, mais même s’il l’était, le résultat avec pivot ne correspondra pas directement au résultat sans pivot. Dans le contexte d’une solution pivotante, ce que Alpha vous a réellement donné est probablement l’ équivalent de ceci:

      1 0 0 1 0 0 1 3 -2 P= 0 1 0 L= 2 1 0 U = 0 -2 7 0 0 1 3 2 1 0 0 -2 

    qui satisfera alors à la condition A = PLU (où. désigne le produit de la masortingce). Si je calcule la même opération de décomposition (théoriquement) d’une autre manière (en utilisant la routine LAPACK dgetrf via numpy dans ce cas):

     In [27]: A Out[27]: array([[ 1, 3, -2], [ 3, 5, 6], [ 2, 4, 3]]) In [28]: import scipy.linalg as la In [29]: LU,ipivot = la.lu_factor(A) In [30]: print LU [[ 3. 5. 6. ] [ 0.33333333 1.33333333 -4. ] [ 0.66666667 0.5 1. ]] In [31]: print ipivot [1 1 2] 

    Après un peu de magie noire avec ipivot nous obtenons

      0 1 0 1 0 0 3 5 6 P = 0 0 1 L = 0.33333 1 0 U = 0 1.3333 -4 1 0 0 0.66667 0.5 1 0 0 1 

    qui satisfait également A = PLU. Ces deux factorisations sont correctes, mais elles sont différentes et ne correspondent pas à une version du code NR qui fonctionne correctement.

    Donc, avant de pouvoir décider si vous avez la “bonne” réponse, vous devriez vraiment passer un peu de temps à comprendre l’algorithme réel implémenté par le code que vous avez copié.