(* file NUMERIAL.M *) (* :Title: NUMERIAL *) (* :Context: MathSource` *) (* :Author: Jesus Rojo *) (* :Summary: Contains general direct methods for solving linear systems *) (* :Mathematica Version: 2.0 *) BeginPackage["NUMERIAL`","NUMERICO`"] ALFactor::usage="ALFactor[A,b,n,prec,lu,pi,pp] decomposes the matrix A (with or without pivoting) in the form L U , and obtains the solution to the system A x = b . The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth \"l\" or \"u\" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth \"cp\" or \"sp\" and indicates if the method is carried out with or without partial pivoting. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The output will be printed in the screen and also in a file called LU.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows and the solution X of the system remain in the memory." ALCholeski::usage="ALCholeski[A,b,n,prec,pp] decomposes the matrix A in the form L tL , and obtains the solution to the system A x = b . The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The output will be printed in the screen and also in a file called CH.DAT. The matrix L , the rounded matrix A (under the name AR) and the solution X of the system remain in the memory." ALFactorwithRefinement::usage= "ALFactorwithRefinement[A,b,n,prec,lu,pi,no,tol,pp,pr] decomposes the matrix A (with or without pivoting) in the form L U , obtains the solution to the system A x = b and improves it through refinement until the relative difference is less than tol. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth \"l\" or \"u\" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth \"cp\" or \"sp\" and indicates if the method is carried out with or without partial pivoting. The real tol gives us the criteria for ending the iterative refinement. The string no is worth \"n1\" , \"n2\" or \"ni\" and indicates that the norm to be compared with tol is either the 1-norm, the euclidean norm or the maximum one. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The string pr is worth \"dr\" or \"sr\" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed in the screen and also in a file called LR.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows, the solution X of the system and the succesive refinements XR remain in the memory." ALCholeskiwithRefinement::usage= "ALCholeskiwithRefinement[A,b,n,prec,no,tol,pp,pr] decomposes the matrix A in the form L tL , obtains the solution to the system A x = b and improves it through refinement until the relative difference is less than tol. The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The real tol gives us the criteria for ending the iterative refinement. The string no is worth \"n1\" , \"n2\" or \"ni\" and indicates that the norm to be compared with tol is either the 1-norm, the euclidean norm or the maximum one. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The string pr is worth \"dr\" or \"sr\" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed in the screen and also in a file called CR.DAT. The matrix L , the rounded matrix A (under the name AR), the solution X of the system and the succesive refinements XR remain in the memory." ALFactorwithnRefinements::usage= "ALFactorwithnRefinements[A,b,n,prec,lu,pi,nref,pp,pr] decomposes the matrix A (with or without pivoting) in the form L U , obtains the solution to the system A x = b and improves it through nref refinements. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The character lu is worth \"l\" or \"u\" and indicates the matrix (L or U) that conserves the diagonal; the other possesses a diagonal equal to 1. The string pi is worth \"cp\" or \"sp\" and indicates if the method is carried out with or without partial pivoting. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The string pr is worth \"dr\" or \"sr\" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed in the screen and also in a file called LN.DAT. The matrices L and U , the matrix P.A or the rounded A (both under the name PA), the vector P of permutation of the rows, the solution X of the system and the succesive refinements XR remain in the memory." ALCholeskiwithnRefinements::usage= "ALCholeskiwithnRefinements[A,b,n,prec,nref,pp,pr] decomposes the matrix A in the form L tL , obtains the solution to the system A x = b and improves it through nref refinements. The method checks if the matrix A is symmetric and if it is positive definite. The integer n is the size of the square matrix A . The integer prec is the number of significant digits with which one wishes to work the method. The string pp is worth \"dp\" or \"sp\" and indicates if the method is carried out with or without partial double precision. The string pr is worth \"dr\" or \"sr\" and indicates if the method is carried out with or without double precision in the refinement. The output will be printed in the screen and also in a file called CN.DAT. The matrix L , the rounded matrix A (under the name AR), the solution X of the system and the succesive refinements XR remain in the memory." Begin["`Private`"] ALFactor[A_,b_,n_Integer,prec_Integer,lu_String,pi_String,pp_String]:= Module[{a,c,B,C,Pivote,pivote,npivote,almacen,divisor,s,i,j,k,fa}, OpenWrite["LU.DAT"]; WriteString["LU.DAT","\n"]; Print[""]; If[pp=="sp"||pp=="SP",fa=1, If[pp=="dp"||pp=="DP",fa=2,Goto[novale]]]; a=0;c=0; If[lu=="u"||lu=="U",Goto[esu],]; If[lu=="l"||lu=="L",Goto[esl],Goto[novale]]; (**) Label[esu]; a=0; Goto[pruebapi]; (**) Label[esl]; a=1; Goto[pruebapi]; (**) Label[pruebapi]; If[pi=="cp"||pi=="CP",Goto[escp],]; If[pi=="sp"||pi=="SP",Goto[essp],Goto[novale]]; (**) Label[escp]; c=0; Goto[general1]; (**) Label[essp]; c=1; Goto[general1]; (**) Label[novale]; WriteString["LU.DAT","\n"]; WriteString["LU.DAT","OPTIONS lu / pi / pp NOT VALID\n"]; WriteString["LU.DAT","\n"]; WriteString["LU.DAT","\n"]; Print[""]; Print["OPTIONS lu / pi / pp NOT VALID"]; Print[""]; Print[""]; Goto[fin]; (**) Label[general1]; B=Table[0.,{n},{n+1}]; C=Table[0,{n},{n+1}]; PA=Table[0.,{n},{n}]; L=Table[0,{n},{n}]; U=Table[0,{n},{n}]; P=Table[i,{i,n}]; X=Table[0.,{n}]; Pivote=Table[0.,{n}]; Do[ Do[ C[[i,j]]=A[[i,j]] ,{i,1,n}] ,{j,1,n}]; Do[ C[[i,n+1]]=b[[i]] ,{i,1,n}]; Do[ Do[ B[[i,j]]=CORounding[A[[i,j]],prec] ,{i,1,n}] ,{j,1,n}]; Do[ B[[i,n+1]]=CORounding[b[[i]],prec] ,{i,1,n}]; WriteString["LU.DAT","For the matrix\n"]; Print["For the matrix"]; WriteString["LU.DAT","\n"]; Print[""]; WriteString["LU.DAT","[ A | b ] = \n"]; Write["LU.DAT",TextForm@MatrixForm[C]]; Print["[ A | b ] = ",MatrixForm[C]]; WriteString["LU.DAT","\n"]; Print[""]; WriteString["LU.DAT","of size ",n," that is\n"]; Print["of size ",n," that is"]; WriteString["LU.DAT","\n"]; Print[""]; WriteString["LU.DAT","[ A | b ] = \n"]; Write["LU.DAT",TextForm@MatrixForm[N[B,prec]]]; Print["[ A | b ] = ",MatrixForm[N[B,prec]]]; WriteString["LU.DAT","\n"]; Print[""]; WriteString["LU.DAT","with ",prec," significant decimal places, the decomposition L U turns out to be\n"]; Print["with ",prec," significant decimal places, the decomposition L U turns out to be"]; WriteString["LU.DAT","\n"]; Print[""]; If[fa==2, WriteString["LU.DAT","(with partial double precision)\n"]; Print["(with partial double precision)"]; WriteString["LU.DAT","\n"]; Print[""] ,]; C=Table[0.,{n},{n}]; Do[ Do[ C[[i,j]]=B[[i,j]] ,{i,1,n}] ,{j,1,n}]; (* son los calculos *) If[c==0,Goto[unoccero],]; If[c==1,Goto[unocuno],]; (**) Label[unoccero]; Do[ Pivote[[i]]=B[[P[[i]],1]] ,{i,1,n}]; pivote=0.; npivote=1; Do[ If[pivotetol, (* ref refinamientos *) m++, almacen=X; Do[ (* X=Pb-PA.almacen; *) s=0.; Do[ s=CORounding[s+CORounding[PA[[i,j]]*almacen[[j]],fb*prec],fb*prec] ,{j,1,n}]; X[[i]]=CORounding[Pb[[i]]-s,fb*prec] ,{i,1,n}]; (* resolver L y = r (o P.r) *) X[[1]]=CORounding[X[[1]]/L[[1,1]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[k,i]]*X[[i]],prec],prec] ,{i,1,k-1}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/L[[k,k]],prec] ,{k,2,n}]; (* fin de resolver L y= r (o P.r) *) (* resolver U x = y *) X[[n]]=CORounding[X[[n]]/U[[n,n]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[U[[k,i]]*X[[i]],prec],prec] ,{i,k+1,n}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/U[[k,k]],prec] ,{k,n-1,1,-1}]; (* fin de resolver U x = y *) Do[ X[[i]]=CORounding[almacen[[i]]+X[[i]],prec] ,{i,1,n}]; XR[[m+1]]=X; If[nno==1,Goto[no1]]; If[nno==2,Goto[no2]]; If[nno==3,Goto[no3]]; (**) Label[no1]; norma=CORounding[CONorm1[X-almacen,n,prec]/CONorm1[X,n,prec],prec]; Goto[general8]; (**) Label[no2]; norma=CORounding[CONorm2[X-almacen,n,prec]/CONorm2[X,n,prec],prec]; Goto[general8]; (**) Label[no3]; norma=CORounding[CONormi[X-almacen,n,prec]/CONormi[X,n,prec],prec]; Goto[general8]; (**) Label[general8]; ]; XR=Table[XR[[i,j]],{i,m},{j,n}]; WriteString["LR.DAT","and the refined solution is\n"]; Print["and the refined solution is"]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","X = ",N[X,prec],"\n"]; Print["X = ",N[X,prec]]; WriteString["LR.DAT","\n"]; Print[""]; If[m==ref+1, WriteString["LR.DAT","(the maximum of ",ref," refinements reached)\n"]; WriteString["LR.DAT","\n"]; Print["(the maximum of ",ref," refinements reached)"]; Print[""] ,]; WriteString["LR.DAT","with refinements\n"]; Print["with refinements"]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","XR = \n"]; Write["LR.DAT",TextForm@MatrixForm[N[XR,prec]]]; Print["XR = ",MatrixForm[N[XR,prec]]]; WriteString["LR.DAT","\n"]; Print[""]; If[fb==2, WriteString["LR.DAT","(the refinement with double precision)\n"]; Print["(the refinement with double precision)"]; WriteString["LR.DAT","\n"]; Print[""] ,]; Goto[fin]; (**) Label[noposiblesinpivote]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","decomposition impossible without pivot element\n"]; Print["decomposition impossible without pivot element"]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","\n"]; Print[""]; Goto[fin]; (**) Label[noposibleconpivote]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","singular matrix (decomposition impossible)\n"]; Print["singular matrix (decomposition impossible)"]; WriteString["LR.DAT","\n"]; Print[""]; WriteString["LR.DAT","\n"]; Print[""]; Goto[fin]; (**) Label[fin]; Global`L=NUMERIAL`Private`L; Global`U=NUMERIAL`Private`U; Global`PA=NUMERIAL`Private`PA; Global`P=NUMERIAL`Private`P; Global`X=NUMERIAL`Private`X; Global`XR=NUMERIAL`Private`XR; Close["LR.DAT"] ] ALCholeskiwithRefinement[A_,b_,n_Integer,prec_Integer,no_String,tol_Real ,pp_String,pr_String]:= Module[{nno,test,C,AA,s,i,j,k,almacen,norma,ref,m,fa,fb}, OpenWrite["CR.DAT"]; WriteString["CR.DAT","\n"]; Print[""]; If[no=="n1"||no=="N1",nno=1, If[no=="n2"||n0=="N2",nno=2, If[no=="ni"||n0=="NI",nno=3,Goto[novale]]]]; If[pp=="sp"||pp=="SP",fa=1, If[pp=="dp"||pp=="DP",fa=2,Goto[novale]]]; If[pr=="sr"||pp=="SR",fb=1, If[pr=="dr"||pr=="DR",fb=2,Goto[novale]]]; C=Table[0,{n},{n+1}]; AA=Table[0.,{n},{n+1}]; AR=Table[0.,{n},{n}]; L=Table[0,{n},{n}]; X=Table[0.,{n}]; Do[ Do[ C[[i,j]]=A[[i,j]] ,{i,1,n}] ,{j,1,n}]; Do[ C[[i,n+1]]=b[[i]] ,{i,1,n}]; Do[ Do[ AA[[i,j]]=CORounding[A[[i,j]],prec] ,{j,1,n}] ,{i,1,n}]; Do[ AA[[i,n+1]]=CORounding[b[[i]],prec] ,{i,1,n}]; Do[ Do[ AR[[i,j]]=CORounding[A[[i,j]],prec] ,{j,1,n}] ,{i,1,n}]; test=0; Do[ Do[ If[AA[[i,j]]!=AA[[j,i]],test=1,] ,{j,1,i-1}] ,{i,1,n}]; If[test==1,Goto[nosimetrica],]; Goto[general1]; (**) Label[novale]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","OPTIONS no / pp / pr NOT VALID\n"]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","\n"]; Print[""]; Print["OPTIONS no / pp / pr NOT VALID"]; Print[""]; Print[""]; Goto[fin]; (**) Label[nosimetrica]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","NON SYMMETRIC MATRIX\n"]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","\n"]; Print[""]; Print["NON SYMMETRIC MATRIX"]; Print[""]; Print[""]; Goto[fin]; (**) Label[nodefinidapositiva]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","NOT POSITIVE DEFINITE MATRIX\n"]; WriteString["CR.DAT","\n"]; WriteString["CR.DAT","\n"]; Print[""]; Print["NOT POSITIVE DEFINITE MATRIX"]; Print[""]; Print[""]; Goto[fin]; (**) Label[general1]; WriteString["CR.DAT","For the matrix\n"]; Print["For the matrix"]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","[ A | b ] = \n"]; Write["CR.DAT",TextForm@MatrixForm[C]]; Print["[ A | b ] = ",MatrixForm[C]]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","of size ",n," that is\n"]; Print["of size ",n," that is"]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","[ A | b ] = \n"]; Write["CR.DAT",TextForm@MatrixForm[N[AA,prec]]]; Print["[ A | b ] = ",MatrixForm[N[AA,prec]]]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","with ",prec," significant decimal places, the matrix L of A = L.tL turns out to be\n"]; Print["with ",prec," significant decimal places, the matrix L of A = L.tL turns out to be"]; WriteString["CR.DAT","\n"]; Print[""]; If[fa==2, WriteString["CR.DAT","(with partial double precision)\n"]; Print["(with partial double precision)"]; WriteString["CR.DAT","\n"]; Print[""] ,]; (* son los calculos *) s=AA[[1,1]]; If[s<=10^(-20),Goto[nodefinidapositiva],]; L[[1,1]]=CORounding[Sqrt[s],prec]; Do[ L[[i,1]]=CORounding[AA[[i,1]]/L[[1,1]],prec] ,{i,2,n}]; nodefinida=0; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[k,j]]*L[[k,j]],fa*prec],fa*prec] ,{j,1,k-1}]; s=CORounding[AA[[k,k]]-s,fa*prec]; If[s<=10^(-20),Goto[finnodefinida],]; L[[k,k]]=CORounding[Sqrt[s],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[i,j]]*L[[k,j]],fa*prec],fa*prec] ,{j,1,k-1}]; L[[i,k]]=CORounding[CORounding[AA[[i,k]]-s,fa*prec]/L[[k,k]],prec] ,{i,k+1,n}]; Goto[finbucle]; Label[finnodefinida]; k=n-1; nodefinida=1; Label[finbucle]; ,{k,2,n-1}]; If[nodefinida==1,Goto[nodefinidapositiva],]; s=0.; Do[ s=CORounding[s+CORounding[L[[n,j]]*L[[n,j]],fa*prec],fa*prec] ,{j,1,n-1}]; s=CORounding[AA[[n,n]]-s,fa*prec]; If[s<=10^(-20),Goto[nodefinidapositiva],]; L[[n,n]]=CORounding[Sqrt[s],prec]; (* fin de los calculos *) WriteString["CR.DAT","L =\n",TextForm@MatrixForm[N[L,prec]],"\n"]; Print["L = ",MatrixForm[N[L,prec]]]; WriteString["CR.DAT","\n"]; Print[""]; Do[ X[[i]]=AA[[i,n+1]] ,{i,1,n}]; (* resolver L y = b *) X[[1]]=CORounding[X[[1]]/L[[1,1]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[k,i]]*X[[i]],prec],prec] ,{i,1,k-1}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/L[[k,k]],prec] ,{k,2,n}]; (* fin de resolver L y = b *) (* resolver tL x = y *) X[[n]]=CORounding[X[[n]]/L[[n,n]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[i,k]]*X[[i]],prec],prec] ,{i,k+1,n}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/L[[k,k]],prec] ,{k,n-1,1,-1}]; (* fin de resolver tL x = y *) almacen=Table[0.,{n}]; XR=Table[0.,{101},{n}]; XR[[1]]=X; norma=Infinity; ref=100; For[ m=1, m=!=ref+1&&norma>tol, (* ref refinamientos *) m++, almacen=X; Do[ (* X=b-AA.almacen; *) s=0.; Do[ s=CORounding[s+CORounding[AA[[i,j]]*almacen[[j]],fb*prec],fb*prec] ,{j,1,n}]; X[[i]]=CORounding[AA[[i,n+1]]-s,fb*prec] ,{i,1,n}]; (* resolver L y = r *) X[[1]]=CORounding[X[[1]]/L[[1,1]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[k,i]]*X[[i]],prec],prec] ,{i,1,k-1}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/L[[k,k]],prec] ,{k,2,n}]; (* fin de resolver L y= r *) (* resolver tL x = y *) X[[n]]=CORounding[X[[n]]/L[[n,n]],prec]; Do[ s=0.; Do[ s=CORounding[s+CORounding[L[[i,k]]*X[[i]],prec],prec] ,{i,k+1,n}]; X[[k]]=CORounding[CORounding[X[[k]]-s,prec]/L[[k,k]],prec] ,{k,n-1,1,-1}]; (* fin de resolver tL x = y *) Do[ X[[i]]=CORounding[almacen[[i]]+X[[i]],prec] ,{i,1,n}]; XR[[m+1]]=X; If[nno==1,Goto[no1]]; If[nno==2,Goto[no2]]; If[nno==3,Goto[no3]]; Label[no1]; norma=CORounding[CONorm1[X-almacen,n,prec]/CONorm1[X,n,prec],prec]; Goto[general8]; Label[no2]; norma=CORounding[CONorm2[X-almacen,n,prec]/CONorm2[X,n,prec],prec]; Goto[general8]; Label[no3]; norma=CORounding[CONormi[X-almacen,n,prec]/CONormi[X,n,prec],prec]; Goto[general8]; Label[general8]; ]; XR=Table[XR[[i,j]],{i,m},{j,n}]; WriteString["CR.DAT","and the refined solution is\n"]; Print["and the refined solution is"]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","X = ",N[X,prec],"\n"]; Print["X = ",N[X,prec]]; WriteString["CR.DAT","\n"]; Print[""]; If[m==ref+1, WriteString["CR.DAT","(the maximum of ",ref," refinements reached)\n"]; WriteString["CR.DAT","\n"]; Print["(the maximum of ",ref," refinements reached)"]; Print[""] ,]; WriteString["CR.DAT","with refinements\n"]; Print["with refinements"]; WriteString["CR.DAT","\n"]; Print[""]; WriteString["CR.DAT","XR = \n"]; Write["CR.DAT",TextForm@MatrixForm[N[XR,prec]]]; Print["XR = ",MatrixForm[N[XR,prec]]]; WriteString["CR.DAT","\n"]; Print[""]; If[fb==2, WriteString["CR.DAT","(the refinement with double precision)\n"]; Print["(the refinement with double precision)"]; WriteString["CR.DAT","\n"]; Print[""] ,]; Goto[fin]; (**) Label[fin]; Global`L=NUMERIAL`Private`L; Global`AR=NUMERIAL`Private`AR; Global`X=NUMERIAL`Private`X; Global`XR=NUMERIAL`Private`XR; Close["CR.DAT"] ] ALFactorwithnRefinements[A_,b_,n_Integer,prec_Integer,lu_String,pi_String ,nref_Integer,pp_String,pr_String]:= Module[{nno,a,c,B,C,Pivote,pivote,npivote,almacen,divisor,s,i,j,k ,Pb,m,fa,fb}, OpenWrite["LN.DAT"]; WriteString["LN.DAT","\n"]; Print[""]; If[pp=="sp"||pp=="SP",fa=1, If[pp=="dp"||pp=="DP",fa=2,Goto[novale]]]; If[pr=="sr"||pp=="SR",fb=1, If[pr=="dr"||pr=="DR",fb=2,Goto[novale]]]; a=0;c=0; If[lu=="u"||lu=="U",Goto[esu],]; If[lu=="l"||lu=="L",Goto[esl],Goto[novale]]; (**) Label[esu]; a=0; Goto[pruebapi]; (**) Label[esl]; a=1; Goto[pruebapi]; (**) Label[pruebapi]; If[pi=="cp"||pi=="CP",Goto[escp],]; If[pi=="sp"||pi=="SP",Goto[essp],Goto[novale]]; (**) Label[escp]; c=0; Goto[general1]; (**) Label[essp]; c=1; Goto[general1]; (**) Label[novale]; WriteString["LN.DAT","\n"]; WriteString["LN.DAT","OPTIONS lu / pi / pp / pr NOT VALID\n"]; WriteString["LN.DAT","\n"]; WriteString["LN.DAT","\n"]; Print[""]; Print["OPTIONS lu / pi / pp / pr NOT VALID"]; Print[""]; Print[""]; Goto[fin]; (**) Label[general1]; B=Table[0.,{n},{n+1}]; C=Table[0,{n},{n+1}]; PA=Table[0.,{n},{n}]; L=Table[0,{n},{n}]; U=Table[0,{n},{n}]; P=Table[i,{i,n}]; X=Table[0.,{n}]; Pivote=Table[0.,{n}]; Do[ Do[ C[[i,j]]=A[[i,j]] ,{i,1,n}] ,{j,1,n}]; Do[ C[[i,n+1]]=b[[i]] ,{i,1,n}]; Do[ Do[ B[[i,j]]=CORounding[A[[i,j]],prec] ,{i,1,n}] ,{j,1,n}]; Do[ B[[i,n+1]]=CORounding[b[[i]],prec] ,{i,1,n}]; WriteString["LN.DAT","For the matrix\n"]; Print["For the matrix"]; WriteString["LN.DAT","\n"]; Print[""]; WriteString["LN.DAT","[ A | b ] = \n"]; Write["LN.DAT",TextForm@MatrixForm[C]]; Print["[ A | b ] = ",MatrixForm[C]]; WriteString["LN.DAT","\n"]; Print[""]; WriteString["LN.DAT","of size ",n," that is\n"]; Print["of size ",n," that is"]; WriteString["LN.DAT","\n"]; Print[""]; WriteString["LN.DAT","[ A | b ] = \n"]; Write["LN.DAT",TextForm@MatrixForm[N[B,prec]]]; Print["[ A | b ] = ",MatrixForm[N[B,prec]]]; WriteString["LN.DAT","\n"]; Print[""]; WriteString["LN.DAT","with ",prec," significant decimal places, the decomposition L U turns out to be\n"]; Print["with ",prec," significant decimal places, the decomposition L U turns out to be"]; WriteString["LN.DAT","\n"]; Print[""]; If[fa==2, WriteString["LN.DAT","(with partial double precision)\n"]; Print["(with partial double precision)"]; WriteString["LN.DAT","\n"]; Print[""] ,]; C=Table[0.,{n},{n}]; Do[ Do[ C[[i,j]]=B[[i,j]] ,{i,1,n}] ,{j,1,n}]; (* son los calculos *) If[c==0,Goto[unoccero],]; If[c==1,Goto[unocuno],]; (**) Label[unoccero]; Do[ Pivote[[i]]=B[[P[[i]],1]] ,{i,1,n}]; pivote=0.; npivote=1; Do[ If[pivote