ΕΠΙΣΤΗΜΟΝΙΚΟΣ ΥΠΟΛΟΓΙΣΜΟΣ Ι

Ονοματεπώνυμο: Αραβανής Κων/νος  Α .Μ.:3628

1η  Εργαστηριακή άσκηση

 

Μέρος 0

Αρχικά θα περιγράψουμε το σύστημα στο οποίο τρέξαμε τα πειράματα μας:

Το σύστημα μας είχε το λειτουργικό Microsoft Windows Vista Business και ο επεξεργαστής

του ήταν ένας Intel Core Duo T2300 στα 1,66GHz. Είχαμε μία ιεραρχεία μνήμης δύο

επιπέδων με το πρώτο επίπεδο να αποτελείτε από δύο κρυφές μνήμες, η μία για data (L1 DCache)

και η άλλη για inst. (L1 D-Cache) και οι δύο των 32ΚBytes και με τα ακόλουθα

χαρακτηριστικά: 8 way set-associative, 64-byte line size. Η κύρια μνήμη (L2 Cache) μας

είχε μέγεθος 2048ΚBytes και ακριβώς τα ίδια χαρακτηριστικά με τις κρυφές. Όσον αφορά το

λογισμικό στο οποίο εκτελέστηκαν τα προγράμματα μας ήταν η Matlab 7.4 Η βιβλιοθήκη τώρα που χρησιμοποιείται από την Matlab για τις πραξεις BLASS συφωνα με το link:

http://www.mathworks.com/support/solutions/data/1-JDIO3.html?solution=1-JDIO3

είναι η Math Kernel Library (MKL) BLAS για τους Intel επεξεργαστές.

 

Μέρος Α

Στο μέρος αυτό της άσκησης μας ζητήθηκε να υλοποιήσουμε τον αλγόριθμο Strassen (με δική μας υλοποίηση ή από το διαδύκτιο). Τον αλγόριθμο τον βρήκα στο internet  στο παρακάτω link:

 

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=2360

 

Ο αλγόριθμος αυτός φαίνεται παρακάτω:

 

function C = strassen_3628(A, B, nmin)

%strassen_3628  Strassen's fast matrix multiplication algorithm.

%               C = strassen_3628(A, B, NMIN), where A and B are matrices of %               dimension

%               a power of 2, computes the product C = A*B.

%               Strassen's algorithm is used recursively until dimension <= % %               NMIN

%               is reached, at which point standard multiplication is used.

%               The default is NMIN = 8 (which minimizes the total number of

%               operations).

 

%               Reference:

%               V. Strassen, Gaussian elimination is not optimal,

%               Numer. Math., 13 (1969), pp. 354-356.

 

if nargin < 3, nmin = 8; end

 

n = length(A);

if n ~= 2^( log2(n) )

   error('The matrix dimension must be a power of 2.')

end

 

if n <= nmin

   C = A*B;

else

   m = n/2; i = 1:m; j = m+1:n;

   P1 = strassen_3628( A(i,i)+A(j,j), B(i,i)+B(j,j), nmin);

   P2 = strassen_3628( A(j,i)+A(j,j), B(i,i), nmin);

   P3 = strassen_3628( A(i,i), B(i,j)-B(j,j), nmin);

   P4 = strassen_3628( A(j,j), B(j,i)-B(i,i), nmin);

   P5 = strassen_3628( A(i,i)+A(i,j), B(j,j), nmin);

   P6 = strassen_3628( A(j,i)-A(i,i), B(i,i)+B(i,j), nmin);

   P7 = strassen_3628( A(i,j)-A(j,j), B(j,i)+B(j,j), nmin);

   C = [ P1+P4-P5+P7  P3+P5;  P2+P4  P1+P3-P2+P6 ];

end

 

 

strassen_3628.m

 

Ο αλγόριθμος αυτός περιέχει τις προϋποθέσεις που ζητούνται στην άσκηση. Συγκεκριμένα επιτρέπει τον πολλαπλασιασμό μητρώων οποιασδήποτε διάστασης παρέχοντας συγχρόνως την δυνατότητα ο χρήστης να ορίζει το ελάχιστο μέγεθος των μητρώων από το οποίο και κάτω ο πολλαπλασιασμός θα γίνεται με τον κλασσικό τρόπο δηλαδή απευθείας στην Matlab, με βάση την τιμή που θα τεθεί το κατώφλι. Επίσης σε περίπτωση που ο χρήστης δεν ορίσει την τιμή κατωφλιού αυτό τίθεται αυτόματα στην τιμή 8.

 

1.Στο ερώτημα αυτό συγκρίνουμε την mtimes με την strassen_3628. Συγκεκριμένα μετράμε το κόστος της  strassen_3628 σε σχεση με το μέσο κόστος υπολογισμού για την mtimes σε μητρώα αντίστοιχων διαστάσεων. Τα μητρώα πάνω στα οποία γίναν τα πειπάματα είναι πυκνά διαστάσεων: 64x64, 128x128, 264x264, 512x512, 1024x1024 . Η κάθε μέτρηση έγινε για τιμές κατωφλιού 8, 16, 32, 64.

 

Ο κώδικας που χρησιμοποιήσαμε είναι ο ακόλουθος:

 

%arxikopoiisi tixaion mitroon C kai D diastaseon 1024x1024

C=rand(1024,1024);

D=rand(1024,1024);

for j=3:1:6

    k=2^j;

    for i=4:-1:0

        %arxikopoiisi mitroon A kai B diastaseon sto 1o loop 64x64

        %                                        sto 2o loop 128x128

        %                                        sto 3o loop 264x264

        %                                        sto 4o loop 512x512

        %                                        sto 5o loop 1024x1024

        A=C(1:2^i:1024,1:2^i:1024);

        B=D(1:2^i:1024,1:2^i:1024);

        %opos kai edo etsi kai sti klisi tis strassen ektelo mia fora prota tin

        %praksi kathos panta i proti ektelesei epistrefei paraplanitika

        %apotelesmata kai meta ksekino tin xronometrisi tis polles fores gia na

        %paro ton meso oro

        mtimes(A,B);

        tic;

        for l=1:1:25

            mtimes(A,B);

        end

        Mtimes=toc/25

        strassen_3628(A,B,k);

        tic;

        for l=1:1:25

            strassen_3628(A,B,k);

        end

        %sto simeio afto ginetai i metrisi tou kostous tis strassen os

        %pollaplasio tou mesou kostous ipologismou gia ti mtimesantistoixon

        %diastaseon

        Strassen_3628=(toc/25)/Mtimes

    end

end

 

merosA_1.m

 

Συγκεντρωτικά τα αποτελέσματα των μετρήσεων μου φαίνονται στον παρακάτω πίνακα (πρόκειτει στην ουσία για τον μέσο χρόνο της strassen δια αυτόν της mtimes):

 

k/διαστάσεις      μητρώων

64x64

128x128

256x256

512x512

1024x1024

8

34.9834

31.2563

27.8671

25.9118

23.3267

16

5.8894

5.9974

5.7941

5.4859

5.1268

32

2.1558

1.9200

2.0952

2.0999

2.0342

64

1.1402

1.2160

1.3560

1.4574

1.4327

merosA_1.txt

 

Απο τον παραπάνω πίνακα πίνακα φαίνεται οτι όσο μικροτερο το κατώφλι τόσο πιο πολύ χρόνο διαρκεί η strassen από την mtimes. Αυτό συμβαίνει γιατί όσο μικρότερο είναι το κατώφλι τόσο περισσότερες φορές εκτελείται η strassen αναδρομικά άρα τόσο περισσότερος χρόνος καταναλώνεται.

 

2.Στο σημείο αυτο της άσκησης θέλουμε να αποδελίξουμε οτι το εμπρός σφάλμα του θεωρητικά ταχέος strassen είναι χειρότερο από την κλασσική υλοποίηση, δηλαδη την mtimes.

 

Ο κώδικας που χρησιμοποιήσαμε για αυτό είναι ο ακόλουθος:

 

%meros A ipoerotima 2

 

%stoixeia ta opoia pollaplasiazoume me ta tuxaia mitroa pou arxikopoioume

%parakato oste na dimiourgisoume mitroa A, B  me eksisou megala stoixeia,

%me eksisou mikra, me iso kai me "antidiametrika" antistoixa

a=[1e10;1e-10;1;1e10];

b=[1e10;1e-10;1;1e-10];

for i=1:1:4

    %arxikopoiisi  mitroon A kai B se arithmitiki diplis akriveias

    A=a(i)*rand(256);

    B=b(i)*rand(256);

    %arikopoiisi mitroon C kai D se moni akriveia

    C=single(A);

    D=single(B);

    %ektelesi pollaplasiasmou strassen me k=2 se arithmitiki monis

    %akriveias

    E_2=single(strassen_3628(C,D,2));

    %ektelesi pollaplasiasmou strassen me k=16 se arithmitiki monis

    %akriveias

    E_16=single(strassen_3628(C,D,16));

    %ektelesi pollaplasiasmou mtimes se arithmitiki monis akriveias

    F=single(mtimes(C,D));

    %ektelesi pollaplasiasmou se arithmitiki dilis akriveias

    G=double(mtimes(double(A),double(B)));

    %ipologismos eN kai eC gia strassen monia akriveias me k=2 kai mtimes

    %diplis akriveias

    eN_strassen_2=norm(E_2-G,2)/((eps('single')/2)*norm(G,2))

    eC_strassen_2=max(max(abs(E_2-G)./((eps('single')/2)*abs(G))))

    %ipologismos eN kai eC gia strassen monia akriveias me k=16 kai mtimes

    %diplis akriveias

    eN_strassen_16=norm(E_16-G,2)/((eps('single')/2)*norm(G,2))

    eC_strassen_16=max(max(abs(E_16-G)./((eps('single')/2)*abs(G))))

    %ipologismos eN kai eC gia mtimes monia akriveias kai mtimes diplis

    %akriveias

    eN_mtimes=norm(F-G,2)/((eps('single')/2)*norm(G,2))

    eC_mtimes=max(max(abs(F-G)./((eps('single')/2)*abs(G))))

end

 

 merosA_2.m

 

Συγκεντωτικά τα αποτελέσματα που λάβαμε ήταν τα παρακάτω:

 

Για strassen μονης ακρίβειας με k=2

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΠΡΩΤΟΥ ΠΙΝΑΚΑ

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΔΕΥΤΕΡΟΥ ΠΙΝΑΚΑ

eN

eC

Πολύ μεγάλες

Πολύ μεγάλες

9.7409

980.5700

Πολύ μικρές

Πολύ μικρές

7.5891

664.2161

Κανονικές

Κανονικές

8.6833

865.8161

Πολύ μεγάλες

Πολύ μικρες

9.5754

851.3267

merosA_2.txt

 

Για strassen μονης ακρίβειας με k=16

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΠΡΩΤΟΥ ΠΙΝΑΚΑ

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΔΕΥΤΕΡΟΥ ΠΙΝΑΚΑ

eN

eC

Πολύ μεγάλες

Πολύ μεγάλες

1.8249

122.3961

Πολύ μικρές

Πολύ μικρές

1.6814

114.4679

Κανονικές

Κανονικές

1.7180

110.5917

Πολύ μεγάλες

Πολύ μικρες

1.6396

113.8317

merosA_2.txt

 

Για mtimes μονης ακρίβειας

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΠΡΩΤΟΥ ΠΙΝΑΚΑ

ΤΙΜΕΣ ΣΤΟΙΧΕΙΩΝ ΔΕΥΤΕΡΟΥ ΠΙΝΑΚΑ

eN

eC

Πολύ μεγάλες

Πολύ μεγάλες

0.0973

2.8324

Πολύ μικρές

Πολύ μικρές

0.0934

3.0715

Κανονικές

Κανονικές

0.0965

2.9746

Πολύ μεγάλες

Πολύ μικρες

0.0977

3.1041

merosA_2.txt

 

Από τους παραπάνω πίνακες και από τις τιμές των eN και eC επαληθεύεται ότι το εμπρός σφάλμα του θεωρητικά ταχέος strassen είναι χειρότερο από την κλασσική υλοποίηση και πιο συγκεκριμένα το σφάλμα είναι μεγαλύτερο όταν το κατώφλι είναι ίσο με 2 από ότι όταν είναι 16. Επίσης παρατηρούμε ότι ότι μητρώα και να επιλέξουμε, τα σφάλματα που προκύπτουν απο την εφαρμεγή της κάθε μεθόδου και λέγοντας μεθόδου εννοούμε strassen με κατώφλι 2, strassen με κατώφλι 16 και mtimes, είναι της ίδιας τάξης μεγέθους.

 

Όσον αφορά την επιλογή ενός μητρώου με πολύ μεγάλες τιμές στοιχείων και ενός με πολύ μικρές αυτή έγινε γιατί εκτελώντας πράξεις ανάμεσα σε τέτοια μητρώα η Matlab στην ουσία απαλοίφει το μητρώο με τις μικρές τιμες και επειδή όταν καλείτε η strassen με πολύ μικρό κατώφλι εκτελείτε αναδρομικά ακόμα περισσότερες φορες παρουσιαζονται πολυ μεγαλα σφάλματα eN και eC για την strassen με k=2. Εδώ αξίζει να σημειώσουμε τον οτι έναχαρακτηριστικό της strassen είναι ότι συνδιάζει σε προσθαφερέσεις τμήματα του κάθε πολλαπλασιστή.

 

Μέρος Β
Στο σημείο αυτό της άσκησης κάνουμε τον σχετικό έλεγχο για το αν ο δείκτης κατάστασης είναι ίσος με το c δηλαδή με το s1/sn σύμφωνα με τις οδηγίες που μας δόθηκαν απο την άσκηση για την κατασκευή μητρώου M(n,c).

 

Ο κώδικας που χρησιμοποιήθηκε ήταν ο ακόλουθος:

 

%elegxos an to cond(A) einai iso me to c me vasi tis proipotheseis pou

%etethisan apo tin ekfonisi tis askisis

 

%o elegxos ginetai gia ta codition pou xrisimopoioume paakato stin ipoloipi

%askisi kathos kai gia to analogo n

n=10;

c=[10;1e2;1e5;1e8];

%arxikopoiisi tou s1 me vasi tin protropi tis askisis

s1=1;

for i=1:1:4

    sn=s1/c(i);

    %arxikopoiisi pinaka diastasewn (n-2)x1 me ta stoixeia tou na

    %kimenontai apo to sn os to s1

    d=sn+(s1-sn).*rand(n-2,1);

    D=diag([s1,d',sn]);

    [U,R]=qr(rand(n));

    [V,R]=qr(rand(n));

    A=U*D*V';

    ConditionOfM=cond(A)

End

merosB_elegxos_condition.m

 

Και πράγματι τα αποτελέσματα που λάβαμε για το condition ήταν οι τιμές που είχαμε δώσει στο c όπως φαίνεται και στο παρακάτω screenshot:

merosB_elegxos_condition.jpg

 

Για το επόμενο ζητούμενο της άσκησης, δηλαδή την εύρεση του εμπρός σχετικού σφάλματος καθώς επίσης και του άνω φράγματος για την επίλυση του συστήματος Αx=b για Α τα οποία καθορίζονται από την εκφώνηση, ο κώδικας μας φαίνεται παρακάτω:

 

%meros B--to 2 stin ousia ipoerotima

 

%arxikopoiisi pinakon stous opoious filassontai ta apotelesmata gia to piso

%kai to empros sfalma kathos kai to ano fragma gia ta mitroa M kai G

%antistoixa

piso_sfalma_M=zeros(4,1);

ano_fragma_M=zeros(4,1);

empros_sfalma_M=zeros(4,1);

piso_sfalma_G=zeros(2,1);

ano_fragma_G=zeros(2,1);

empros_sfalma_G=zeros(2,1);

%dimiourgia pinaka e nx1 pou periexei monades os stoixeia

e=ones(n,1);

%arxikopoiisi n me vasi tin ekfonisi kathos kai dimiourgia pinaka c pou

%periexei ola ta c pou ha xrisimopoiithoun stin ektelesi

n=10;

c=[10;1e2;1e5;1e8];

%arxikopoiisi s1 me vasi tin protropi tis askisis

s1=1;

%arxikopoiisi metriti j

j=1;

for i=1:1:4

    %arxikopoiisi sn oste na isxiei oti c=s1/sn

    sn=s1/c(i);

    %arxikopoiisi pinaka diastasewn (n-2)x1 me ta stoixeia tou na

    %kimenontai apo to sn os to s1

    d=sn+(s1-sn).*rand(n-2,1);

    D=diag([s1,d',sn]);

    [U,R]=qr(rand(n));

    [V,R]=qr(rand(n));

    A=U*D*V';

    b=A*e;

    %euresi tou eponomazomenou stin ekfonisi x^

    x=A\b;

    %ipologismos piso kai empros sfalmato kathos kai no fragmatos gia ton

    %pinaka M(n,c) tis ekfonisis

    piso_sfalma_M(i)=norm(b-A*x,2)/(norm(A,2)*norm(x,2)+norm(b,2));

    empros_sfalma_M(i)=norm(e-x)/norm(e);

    ano_fragma_M(i)=2*cond(A)*piso_sfalma_M(i);

    if((c(i)==10)|(c(i)==1e5))

        %dimiourgia tou pinaka G

        G=gfpp(triu(A(1:n-1,1:n-1)));

        b=G*e;

        A=G;

        %euresi tou eponomazomenou stin ekfonisi x^

        x=A\b;

        %ipologismos piso kai empros sfalmato kathos kai no fragmatos gia ton

        %pinaka G tis ekfonisis

        piso_sfalma_G(j)=norm(b-A*x,2)/(norm(A,2)*norm(x,2)+norm(b,2));

        empros_sfalma_G(j)=norm(e-x)/norm(e);

        ano_fragma_G(j)=2*cond(A)*piso_sfalma_G(j);

        j=j+1;

    end

end

%ektiposi empros sfalmatos kai ano_fragmatos gia ton kathe M(n,c) kai ton kathe G

empros_sfalma_M

ano_fragma_M

empros_sfalma_G

ano_fragma_G

merosB.m

 

Τα αποτελέσματα που λάβαμε φείνονται στον παρακάτω πίνακα:

 

Εμπρός σφάλμα

(x 1.0e-010)

Άνω φράγμα

(x 1.0e-009)

M(10,10)

0.000

0.00

M(10,1e2)

0.000

0.00

M(10,1e5)

0.027

0.01

M(10,1e8)

8.380

11.14

G για c=10

0.0003

0.0004

G για c=1e5

0.2086

0.1421

merosB.txt

Μέρος Γ

Στο σημείο αυτό της άσκησης έχουμε σαν σκοπό την εύρεση των Mflop/s για την επίλυση του συστήματος Αx=b με μονη και διπλή ακρίβεια, χρησιμοποιώτας τυχαία πυκνά μητρώα μεγεθών 1000:200:2000 και ακολούθως να τα συγκρίνουμε.

 

Ο κώδικας που χρησιμοποιήθηκε για τον υπολογισμό των Mflop/s είναι αυτός:

%merosC euresi Mflop/s

 

%arxikopoiisi sto 0 pinaka ston opoio tha apothikeutoun ta ipologismena

%Mflop/s

mflops=zeros(6,2);

%arxikopoiisi metriti j

j=1;

for i=1000:200:2000

    %dimiourgia tixaion mitroon diastaseon ixi kai ix1 antistoixa

    A=rand(i);

    b=rand(i,1);

    %ektelesi arxika mias foras tis prakseis x=A\b prin ksekinisoume tin

    %xronmetrisi kathos i proti ektelesi mporei na dosi paraplanitika

    %apotelesmata gia ton xrono kai etsi den tin simperilamvanoume sto meso

    %oro

    x=A\b;

    %euresi mesou xronou ektelesis tis prakseis x=A\b se arithmitiki

    %diplis akriveias

    tic;

    for k=1:1:20

        x=A\b;

    end

    dipli_akriveia=toc/20;

    %arxikopoiisi mitroon C kai d monis akriveias

    C=single(A);

    d=single(b);

    %ektelesi mias foras tis prakseis x=C\d gia ton idio opos kai

    %proigoumenos logo

    x=C\d;

    %euresi mesou xronou ektelesis tis prakseis x=C\d se arithmitiki

    %monis akriveias

    tic;

    for k=1:1:20

        x=single(C\d);

    end

    moni_akriveia=toc/20;

    %ipologismos ton Mflop/s

    mflops(j,1)=((2/3)*i^3+2*i^2)/dipli_akriveia;

    mflops(j,2)=((2/3)*i^3+2*i^2)/moni_akriveia;

    j=j+1;

end

x=[1000;1200;1400;1600;1800;2000];

%ektiposi grafikis parastasis me stoixeia tin ekseliksi ton Mflop/s

%gia prakseis diplis kai monis akriveias apo mitroa diastasis

%1000 x 1000 se mitroa 2000 x 2000

plot(x, mflops(:,1), '.-');

% kanoume hold gia na emfanistoun se koini grafiki parastasi

hold on;

plot(x, mflops(:,2), 'g.-');

legend('dipli akriveia', 'moni akriveia')

hold;

MerosC_1.m

 

Η γραφική παράσταση που λάβαμε σαν αποτέλεσμα ήταν η εξής:

 

merosC_1.fig

 

Όπως ήταν αναμενόμενο άλλωστε η αναπαράσταση μονής ακρίβειας εκτελεί πολύ περρισότερες πράξεις το δευτερόλεπτο σε σχέση με την αναπαράσταση διπλής ακρίβειας και αυτό λόγο της μικρότερης πολυπλοκότητας των πράξεων αφου επεξεργάζονται μικρότερης ακρίβειας δεδομένα. Εύκολα διαπυστώνουμε από την γραφική παράσταση ότι όταν μεγαλώνει και το μέγεθος των μητρώων με τα οποία ασχολούμαστε τοσο μεγαλώνουν, ελαφρώς βέβαια και οι αντίστοιχες τιμές των Mflop/s.

 

Τώρα σαν στόχο έχουμε να υλοποιήσουμε κώδικα στον οποίο να παρατηρούμε τις διαφορές στους χρόνους επίλυσης του συστήματος  Ax=b, με αλγόριθμο επαναλυπτικής εκλέπτυνσης σε μονή ακρίβεια και με την κλασσική υλοποίηση, δηλαδή x=A^(-1)*b,  καθώς επίσης και τις μεταξύ διαφορές των αποτελεσμάτων (των x δηλαδή) με είσοδο διάφορων διαστάσεων τυχαία μητρώα.

 

Ο κώδικας που υλοποιήσαμε φαίνεται παρκάτω:

 

%meros C -- ipoerotima 2

 

%arxikopoiisi pinakon pou apothikeuoun ta apotelesmata gia to poses fores

%einai taxiteri i epanalitiki ekleptinsi apo to mtimes, tis diastaseis ton

%pinakon pou xrisimopoioume, tin maximum diafora gia ton x pou

%ipologizontai me tis dio diaforetikes methodoud, tin minimum, ton meso

%xrono tis epanaliptikis ekleptinsisgia kathe diastasi kai ton meso xrono

%tis mtimes antistoixa

poses_fores_taxtiteros=zeros(10,1);

diastaseis=zeros(10,1);

maximum=zeros(10,1);

minimum=zeros(10,1);

epanaliptiki_ekleptinsi=zeros(10,1);

mtimes=zeros(10,1);

%arxikopoiisi metriti j

j=0;

for n=100:100:1000

    j=j+1;

    %dimiourgia tixaion mitroon A,b

    A=rand(n);

    b=rand(n,1);

    tic;

    for i=1:1:20

        %Algorithmos epanaliptikis ekleptinsis

       

        %diaspasi tou A se mitroa L, U kai P

        [L,U,P]=lu(A);

        %euresi ton L, U kai P se moni akriveia

        U_single=single(U);

        L_single=single(L);

        P_single=single(P);

        %ipologismos tou x se moni akriveia

        x=single(U_single\(L_single\(P_single*single(b))));

        k=1;

        while(k==1||norm(z)<=eps('single')*norm(x))

            r=b-A*x;

            z=single(U_single\(L_single\(P_single*single(r))));

            x=x+z;

            k=k+1;

        end

    end

    %apothikeusi mesou xronou ipologismou tou x me tin methodo tis

    %epanaliptikis ekleptinsis

    epanaliptiki_ekleptinsi(j)=toc/20;

    %ipologismos tou x me ton klassiko pollaplasiasmo tis Matlab, edo to x

    %to onomazo y

    tic;

    for i=1:1:20

        y=inv(A)*b;

    end

    %apothikeusi mesou xronou ipologismou tou x me tin methodo tis mtimes

    mtimes(j)=toc/20;

    %ipologismos diaforas metaksi ton apotelesmaton ton dio methodon

    diafora_metaksi_ton_2_x=y-x;

    %ipologismos tis apolitis timis ton diaforon

    l=abs(diafora_metaksi_ton_2_x);

    %ipologismos maximum fiaforas

    maximum(j)=max(l);

    %ipologismos minimum diaforas

    minimum(j)=min(l);

    %topothetisi ton diastaseon ston pinaka me tis diastaseis

    diastaseis(j)=n;

    %ipologismos tou pilikou mtimes/epanaliptiki_ekleptinsi, oste na

    %diapistosoume poses fores einai taxiteros o algorithmos epanaliptikis

    %ekleptinsis apo tin euresi tou x me tin mtimes

    poses_fores_taxtiteros(j)=mtimes(j)/epanaliptiki_ekleptinsi(j);

end

plot(diastaseis,poses_fores_taxtiteros,'.-')

legend('poses fores taxiteros o Algorithmos Epanaliptikis Ekleptinsis apo afto tis Mtimes')

epanaliptiki_ekleptinsi

mtimes

maximum

minimum

poses_fores_taxtiteros

merosC_2.m

 

Στο παρακάτω σχήμα φαίνεται η γραφική παράσταση του πηλικου του μέσου χρόνου εκτέλεσης της επίλυσης του συστήματος με την κλασσική υλοποιηση προς αυτόν με τον αλγόριθμο επαναληπτικής  εκλέπτυνσης:

 

merosC_2.fig

 

Όπως εύκολα μπορούμε να διαπυστώσουμε με μονή ακρίβεια η λύση επέρχεται πολύ πιο γρήγορα, αναμενόμενο άλλωστε λόγο της πολυπλοκότητας των υπολογισμών σε σχέση με την διπλή ακρίβεια. Επίσης όσο μεγαλύτερα τα μητρώα με τα οποία εχουμε να κάνουμε τόσο μεγαλύτερος είναι και ο λόγος των δύο χρόνων, αν και όσο τίνουμε προς πολυ μεγάλα μητρώα ο λόγος αυτός υποθέτουμε ότι συγκλίνει.

 

Συγκεκριμένα οι αντίστιχοι χρόνοι για την μέθοδο της επαναληπτικής εκλέπτυνσης και της mtimes,από τους οποίους προκύπτει ο παραπάνω λόγος, φαίνονται στον παρακάτω πίνακα:

 

ΔΙΑΣΤΑΣΕΙΣ Α

ΧΡΟΝΟΣ ΓΙΑ ΕΠΑΝΑΛΗΠΤΙΚΗΣ ΕΚΛΕΠΤΥΝΣΗΣ

ΧΡΟΝΟΣ ΓΙΑ MTIMES

100x100

0.0022

0.0028

200x200

0.0071

0.0148

300x300

0.0219

0.0490

400x400

0.0492

0.1119

500x500

0.0902

0.2141

600x600

0.1497

0.3687

700x700

0.2308

0.5779

800x800

0.3394

0.8613

900x900

0.4704

1.2119

1000x1000

0.6347

1.6544

merosC_2.txt

 

Τώρα βέβαια γεννάται η απορία για την απώλια ακριβειας της λύσης x, όταν οι υπολογισμοί γίνονται σε αριθμητική μονής ακρίβειας σε σχεση με αυτή της διπλής. Παρ’ όλα αυτα από τον παρακάτω πίνακα είναι ολοφάνερο οτι η μέγιστη διαφορά στα μεταξύ διπλής και μονής ακρίβειας αποτελέσματα, δηλαδή τα x, είναι αμελιτέα σε σχέση με τον χρόνο που γλιτώνουμε. Τώρα για την ελάχιστη διαφορά δεν αξίζει να γίνει λόγος αφου είναι υπερβολικά μικρή.

 

ΔΙΑΣΤΑΣΕΙΣ Α

MAXIMUM ΔΙΑΦΟΡΑ

MINIMUM ΔΙΑΦΟΡΑ

(x1.0e-005)

100x100

0.0000

0.0059

200x200

0.0000

0.0039

300x300

0.0018

0.0166

400x400

0.0000

0.0008

500x500

0.0000

0.0002

600x600

0.0001

0.0011

700x700

0.0008

0.1037

800x800

0.0001

0.0050

900x900

0.0000

0.0010

1000x1000

0.0004

0.0011

merosC_2.txt

 

Η διαφορές αυτές είναι ανεξάρτητες από το μέγεθος του πίνακα που χρησιμοποιούμε σαν δεδομένο για την επίλυση των συστημάτων. Έπίσης είναι πολύ μικρές λόγω του ότι δεν λύθηκε απλά το σύστημα σε αριθμητική μονής ακρίβειας αλλά υλοποιήθηε ένας αλγοριθμος επαναληπτικής εκλέπτυνσης.