บทที่ 4 สปาร์สเมทริกซ์ และตัวแก้ปัญหาที่ตรง

บทที่ 4

สปาร์สเมทริกซ์ และตัวแก้ปัญหาที่ตรง (Sparse Matrices and Direct Solvers)

4.1 ส่วนนำ

            เอกสารนี้ออกแบบเพื่อสร้างความคุ้นเคยให้ผู้ใช้กับการคำนวณสปาร์สเมทริกซ์ในไซแลบ

4.2 สปาร์สเมทริกซ์

            ไซแลบยอมให้ผู้ใช้ทำงานกับสปาร์สเมทริกซ์ทั้งหมดได้ง่ายๆเท่ากับเมทริกซ์เต็มทั้งหมด ให้มาสร้างสปาร์สเมทริกซบางอย่าง

4.3 การแปลเมทริกซ์ทั้งหมดไปเป็นเมทริกซ์สปาร์ส

            สามารถที่จะใช้คำสั่งdiag เพื่อที่จะสร้างเมทริกซ์ tridiagonal  ให้ทดลอง

A = diag(ones(20,1),-1) -2*diag(ones(21,1)) + diag(ones(20,1),1);

            จากบรรทัดคำสั่งก่อให้เกิดเมทริกซ์ทั้งหมดขนาด 21 × 21  โดยผู้ใช้สามารถแปลงไปสู่พื้นที่เก็บ สปาร์สเมทริกซ์  โดยใช้คำสั่ง sparse ดังนี้

A = sparse(A);

            ด้วยเมทริกซ์สปาร์ส  จะมีประโยชน์ที่จะมองไปที่โครงสร้างของเมทริกซ์ที่ไม่เป็นศูนย์ (non-zero)

หนังสือ“Engineering and Scientific Computing with Scilab” ที่ปรับปรุงโดย Claude Gomez ได้จัดส่วนแนะนำที่มีประโยชน์มากในการใช้งานไซแลบ โดยเฉพาะอย่างยิ่งที่ทำให้มีบางรหัสคำสั่งซึ่งลอกเลียนคำสั่ง spy จาก Matlab. ต่อไปได้สร้างรหัสคำสั่งที่ฝั่งอยู่ในฟังก์ชันดังนี้

//-------------------------------------------------------

function spy(A)

// Mimic the Matlab spy command

// Draw the non zero components of a matrix

// Taken from Engineering and Scientific Computing with Scilab’’

// edited by Claude Gomez

[i,j] = find(A~=0)

[N,M] = size(A)

xsetech([0,0,1,1],[1,0,M+1,N])

xrects([j;N-i+1;ones(i);ones(i)],ones(i));

xrect(1,N,M,N);

endfunction

//-------------------------------------------------------

ให้สำเนารหัสคำสั่งข้างบนเข้าไปในพื้นที่ทำงานไซแลบ และสั่งให้ทำงานกับเมทริกซ์ A

xbasc();spy(A)

            คำสั่ง xbasc()  เป็นเพียงเพื่อล้างวินโดว์กราฟิกส์  ถ้าคุณได้ใช้ไปเรียบร้อยแล้ว

4.4 การสร้างสปาร์ส เมทริกซ์

            เมื่อสร้าง A ในครั้งแรกให้เป็นเมทริกซ์เต็มทั้งหมด แล้วแปลงไปเป็นสปาร์สเมทริกซ์อีกที่  เป็นที่ชัดแจ้งง่าควรจะหลีกเลี่ยงเมทริกซ์ที่เต็มทั้งหมด  ให้ใช้ whos() หรือwhos -name A  เพื่อดูว่าต้องการที่เก็บเท่าใดของ A ก่อนและหลังการแปลงไปอยู่ในรูปแบบการเก็บสปาร์สเมทริกซ์

ตอนนี้ให้ทดลองสร้างเมทริกซ์เป็นแบบสปาร์สเมทริกซ์โดยไม่จำเป็นต้องใช้ส่วนเก็บแบบเต็มทั้งหมด  สามารถที่ใช้ฟังก์ชัน diag อีกครั้งแต่มีอาร์กิวเมนต์หนึ่งเป็นสปาร์ส  คำสั่งต่อไปนี้จะทำกลไกนี้

n=21; B = diag(sparse(ones(n-1,1)),-1) ...

- 2*diag(sparse(ones(n,1))) + diag(sparse(ones(n-1,1)),1);

            ตรวจสอบเมทริกซ์โดยพิมพ์ส่วนที่เป็นข้อมูลนำเข้า หรือแม้แต่พิมพ์แบบเต็มที่แทน B .

ไม่ทำเช่นนี้กับเมทริกซ์ขนาดใหญ่  ให้ใช้ spy ในการตรวจสอบโครงสร้าง รูปแบบการเก็บของสปาร์สโดยหลักๆ แล้วเป็นลิสท์หนึ่งของการนำเข้าองค์ประกอบ ij ที่ตรงกับเมทริกซ์ค่า v ของข้อมูลนำเข้าในเมทริกซ์ที่ไม่เป็นศูนย์ ที่ขณะใดๆ ให้ทดลองทำต่อไปนี้

//-------------------------------------------------------

dij = [(1:n)’, (1:n)’];                  // i j coordinates of diagonal

dv = [-2*ones(n,1)];                  // value of -2 down diagonal

udij = [ (1:n-1)’,(2:n)’];                         // i j coordinates of upper diagonal

udv = [ ones(n-1,1) ];               // value of 1 for upper and lower diagonals

ldij = [ (2:n)’,(1:n-1)’];              // i j coordinates of lower diagonal

ij = [dij ; udij ; ldij ];                 // concatenate all the i, j coordinates

v = [dv ; udv ; udv ];                // concatenate all the values

C=sparse(ij,v)                           // create sparse matrix

full(C)

//-------------------------------------------------------

            ให้แน่ใจว่ารู้ถึงรหัสคำสั่งที่กำลังกระทำ  จงทดลองค่า n ที่น้อยกว่าและพิมพ์ขั้นตอนแต่ละขั้น และให้เปรียบเทียบ เมทริกซ์ A, B, C. จงอธิบายแต่ละแบบสร้างขึ้นอย่างไร ให้แสดงความคิดเห็นถึงประสิทธิภาพของแต่ละวิธี (อัตราเร็วและที่เก็บข้อมูล)

4.5 ตัวอย่างเมทริกซ์ กราฟ

            ไฟล์ความช่วยเหลือสำหรับฟังก์ชัน mesh2d จะได้เมทริกซ์ที่น่าสนใจมากเกี่ยวข้องกับบริเวณพื้นที่สามเหลี่ยยมและวงกลม  ต่อไปจะได้สืบทอดตัวอย่างที่กำหนดไว้ในไฟล์ความช่วยเหลือ  ให้สำเนารหัสคำสั่งนี้เข้าไปสู่พื้นที่ทำงานไซแลบ.

//-------------------------------------------------------

function [a1,a2,a3,g1,g2,g3]=mesh_examples()

//

// Code taken from the help file for mesh2d. Produces 3

// matrices and graphs associated with triangulating

// a circular region with a hole

//

// FIRST CASE

theta=0.025*[1:40]*2.*%pi;

x=1+cos(theta);

y=1.+sin(theta);

theta=0.05*[1:20]*2.*%pi;

x1=1.3+0.4*cos(theta);

Computational Science, Scilab Tutorials 31

y1=1.+0.4*sin(theta);

theta=0.1*[1:10]*2.*%pi;

x2=0.5+0.2*cos(theta);

y2=1.+0.2*sin(theta);

x=[x x1 x2];

y=[y y1 y2];

[nu,a1]=mesh2d(x,y);

nbt=size(nu,2);

jj=[nu(1,:)’ nu(2,:)’;nu(2,:)’ nu(3,:)’;nu(3,:)’ nu(1,:)’];

as=sparse(jj,ones(size(jj,1),1));

ast=tril(as+abs(as’-as));

[jj,v,mn]=spget(ast);

n=size(x,2);

g=make_graph(’foo’,0,n,jj(:,1)’,jj(:,2)’);

g(’node_x’)=300*x;

g(’node_y’)=300*y;

g(’default_node_diam’)=10;

g1 = g;

// SECOND CASE !!! NEEDS x,y FROM FIRST CASE

x3=2.*rand(1:200);

y3=2.*rand(1:200);

wai=((x3-1).*(x3-1)+(y3-1).*(y3-1));

ii=find(wai >= .94);

x3(ii)=[];y3(ii)=[];

wai=((x3-0.5).*(x3-0.5)+(y3-1).*(y3-1));

ii=find(wai <= 0.055);

x3(ii)=[];y3(ii)=[];

wai=((x3-1.3).*(x3-1.3)+(y3-1).*(y3-1));

ii=find(wai <= 0.21);

x3(ii)=[];y3(ii)=[];

xnew=[x x3];ynew=[y y3];

fr1=[[1:40] 1];fr2=[[41:60] 41];fr2=fr2($:-1:1);

fr3=[[61:70] 61];fr3=fr3($:-1:1);

front=[fr1 fr2 fr3];

[nu,a2]=mesh2d(xnew,ynew,front);

nbt=size(nu,2);

jj=[nu(1,:) ’ nu(2,:)’;nu(2,:)’ nu(3,:)’;nu(3,:)’ nu(1,:)’];

as=sparse(jj,ones(size(jj,1),1));
ast=tril(as+abs(as’-as));
[jj,v,mn]=spget(ast);
n=size(xnew,2);
g=make_graph(’foo’,0,n,jj(:,1)’,jj(:,2)’);
g(’node_x’)=300*xnew;
g(’node_y’)=300*ynew;
g(’default_node_diam’)=10;
g2=g;
// REGULAR CASE !!! NEEDS PREVIOUS CASES FOR x,y,front
xx=0.1*[1:20];
yy=xx.*.ones(1,20);
zz= ones(1,20).*.xx;
x3=yy;y3=zz;
wai=((x3-1).*(x3-1)+(y3-1).*(y3-1));
ii=find(wai >= .94);
x3(ii)=[];y3(ii)=[];
wai=((x3-0.5).*(x3-0.5)+(y3-1).*(y3-1));
ii=find(wai <= 0.055);
x3(ii)=[];y3(ii)=[];
wai=((x3-1.3).*(x3-1.3)+(y3-1).*(y3-1));
ii=find(wai <= 0.21);
x3(ii)=[];y3(ii)=[];

xnew=[x x3];ynew=[y y3];

[nu,a3]=mesh2d(xnew,ynew,front);

nbt=size(nu,2);

jj=[nu(1,:)’ nu(2,:)’;nu(2,:)’ nu(3,:)’;nu(3,:)’ nu(1,:)’];

as=sparse(jj,ones(size(jj,1),1));

ast=tril(as+abs(as’-as));

[jj,v,mn]=spget(ast);

n=size(xnew,2);

g=make_graph(’foo’,0,n,jj(:,1)’,jj(:,2)’);

g.node_x=300*xnew;

g(’node_y’)=300*ynew;

g(’default_node_diam’)=3;

g3=g;

endfunction

//--------------------------------------

// Now run the function to produce the

// example matrices and associated graphs

//--------------------------------------

[A1,A2,A3,g1,g2,g3] = mesh_examples();

//-------------------------------------------------------

            หากได้โหลดรหัสคำสั่งที่แล้วเข้ามาในไซแลบ ก็ควรจะมีเมทริกซ์ใหม่ 3 เมทริกซ์ A1, A2, A3 กำหนดไว้ในพื้นที่ทำงานผู้ใช้  ให้ทดลองใช้คำสั่ง spy บนเมทริกซ์เหล่านี้  การได้มาซึ่งการพล็อต spy ของเมทริกซ์ A3  ให้พิมพ์แสดงผลการพล็อต spy ของเมทริกซ์ A2 . คิดว่าน่าจะสร้างความเข้าใจกับการพล็อต spy
            จากนี้รหัสคำสั่งที่แล้วก่อให้เกิดการสร้าง 3 ออบเจ็คของกราฟกับแต่ละสปาร์สเมทริกซ์  จงดูที่กราฟโดยใช้คำสั่งแสดงกราฟ (ถ้าผู้ใช้ต้องการทำแบบนี้บนระบบปฏิบัติการวินโดว์ และใช้ไซแลบรุ่นแรกๆ เช่นรุ่น 2.6 ของไซแลบ จำเป็นต้องโหลดแพคเกจซ์ scigraph เข้ามาก่อนและใช้คำสั่ง show_scigraph แทน)

xbasc(); show_graph(g1);

ปกติแล้วคิดให้ g2 และ g3 เป็นที่น่าสนใจมากกว่า  สามารถใช้ทางเลือกของเมนูย่อยกราฟเพื่อให้จำนวนของโหนด  เรื่องสามารถช่วยให้เข้าใจโครงสร้างที่แปลงของสปาร์สเมทริกซ์ที่ตรงกัน  การได้มาซึงการพล็อตกราฟต่อไปนี้เกี่ยวข้องกับเมทริกซ์ A3 และให้พิมพ์การพล็อตกราฟที่ตรงกับเมทริกซ์ A2 ให้มองไปที่คอลัมน์แรกของเมทริกซ A2.  โดยทำให้ลงรอยกับกราฟ g2

4.6 การดำเนินการเมทริกซ์สปาร์ส

            การดำเนินการเมทริกซ์มาตรฐานเกือบทั้งหมดนำไปใช้กับเมทริกซ์สปาร์ส โดยสามารถ บวก คูณ ทรานสฟอร์ม(transform) อินเวิร์ท(invert) และแก้สมการเมทริกซ์  เช่น สามารถจะหาอินเวิรทของเมทริกซ์ A1.ที่สร้างจากตัวอย่างฟังก์ชัน mesh_examples

A1inv = inv(A1);

            มองไปที่ A1inv กับฟังก์ชัน spy . การอินเวิร์สเกือบจะเต็ม เป็นปกติที่อินเวิร์สของเมทริกซ์สปาร์สแบบสุ่มที่เปลี่ยนไปเกือบเต็ม สำหรับ moral, หลีกเลี่ยงการคำนวณอินเวิร์สของเมทริกซ์สปาร์ส

4.7 การแก้สมการเมทริกซ์

            โดยพื้นฐานแล้วมี 3 วิธีที่หาได้ในการแก้สมการเมทริกซ์

•  ตัวกระทำแบคสแลส (backslash operator \)

• การควบรวม lufact และ lusolve

• และสำหรับเมทริกซ์หาค่าบวกได้ จากการควบรวมฟังก์ชัน chfact และ chsolve

• ฟังก์ชันspchol นำมาใช้ได้แต่ฟังก์ชัน chfact ค่อนข้างมีประสิทธิภาพดีกว่า

ติดตั้งปัญหาทดสอบ ดังนี้

[n,m]=size(A1);

x = ones(n,1);

b = A1*x;

            ดังนั้นมี b ทางด้านขวาร่วมด้วยกับเว็คเตอร์คำตอบ x   จากนี้ดูว่าแต่ละวิธีการที่ได้คำตอบทำงานดีอย่างไร ให้ใช้คำสั่ง help เพื่อดูวิธีการเรียกใช้ factorization ที่แตกต่างกัน  และใช้คำสั่ง timer() เพื่อจับเวลาการคำนวณแต่ละครั้ง  ดังเช่น

timer(); x1 = A1\b ; t1 =timer()

err1 = norm(x1-x)

4.8 ส่วนขยาย Extensions

            คำสั่งควบรวมเช่น chfact, chsolve  ค่อนข้างมีประสิทธิภาพสำหรับเมทริกซ์หาค่าได้เป็นบวกตามที่ใช้องศาการจัดลำดัใหม่ต่ำสุด (สำหรับ lufact และ spcholไม่ได้จัดลำดับใหม่) น่าเสียดายว่าเมทริกซ์ A1 ไม่ได้เป็นเมทริกซ์ที่หาค่าได้เป็นบวก  ติดตามดูว่าสามารถที่จะปรับปรุงประสิทธิภาพของคำสั่ง lufact โดยจัดลำดับเมทริกซ์ A1 ใหม่  คำสั่ง bandwr จะพยายามที่จะจัดลำดับเมทริกซ์สมมาตรใหม่ เพื่อสร้างเมทริกซ์ที่จัดลำดับใหม่ด้วยช่วงแบนวิดท์ที่เล็กกว่า  ประยุกต์วิธีการนี้กับเมทริกซ์ A1 (ซึ่งเป็นเมทริกซ์สมมาตร)  ทดลองดังนี้

[iperm,mrepi]=bandwr(triu(A1));

RA1 = A1(mrepi,mrepi);

xbasc();spy(A1)

xbasc();spy(RA1)

            ข้อสังเกตก็คือการจัดลำดับเมทริกซ์ใหม่มีแบนวิดท์ที่เล็กกว่า ให้หาว่าวิธีการหาคำตอบทำงานได้ดีกว่าหรือไม่ในการจัดลำดับเมทริกซ์ใหม่