Description

Cabi is a student at NTHU, who recently has faced some trouble with linear algebra.

Recently, he was asked to solve for the pressure of a fluid system using the Navier-Stokes equations, which consist of the Mass Conservation equation: and the Momentum Conservation equation: .

Cabi was given a mesh file which described the fluid system, and he found out that he can transform the problem into a system of linear equations. Given a matrix and a vector , his goal is to solve for the unknown vector in .

The values in each row of matrix represent boundary information between neighbors. Since we are in a 3D world, the number of neighbors for each cell is far less than the total number of cells in the system, making the matrix here a sparse matrix.

Due to the nature of the physical properties, the matrix here is a symmetric matrix, which is clearly true.

• A sparse matrix means a matrix with a lot of zeros.
• A symmetric matrix means the coefficient at is equal to the one at .

Cabi had found the matrix , but he ran out of his paper._..

Having no paper and no coding skills, Cabi suffered a mental breakdown.

An idea crossed Cabi's mind. As an I2P TA without empathy, he decided to let every I2P student solve this problem for him, he thought that will be a piece of cake for these geniuses.

---

Now, your task is to help Cabi solve the problem. Write a program to solve the linear system using the Conjugate Gradient(CG) algorithm.

Unlike Gaussian Elimination,CG is an iterative method perfect for large, sparse, and symmetric positive-definite matrices. Here is how it works:

0. Prepare a starting vector called (demension=n), the matrix () and the vector (demension=n).
1. Initiate variables: let
2. Iterate the below for

Note:
In linear algebra, a vector will be a column vector, and means change the row and column, so means .

As a kind TA, Cabi has already written a basic structure for you, but some of the function are still missing. Your task is to complete those functions.

Variables

• N: the number of cells in the fluid system. This is also the size of the vectors and the matrix
• step: the number of iterations you need to perform for the CG algorithm.
• d: a list of length . d[i] represents the diagonal element of matrix at row (i.e.,).
• pos: a list of lists. pos[i] contains the column indices of the non-zero off-diagonal(非對角線) elements in row of the matrix .
• val: a list of lists. val[i] contains the values of the non-zero off-diagonal elements in row , corresponding to pos[i] of the matrix .
• b: the right-hand side of the equation of length .
• x0: the initial guess vector of length .

Example for the sparse matrix format:

Suppose we have a matrix :

The inputs for this matrix will look like this:

4 4 4 (this is d)
1 1 -1           (Row 0: 1 off-diagonal element at column 1, value is -1)
2 0 2 -1 -2      (Row 1: 2 off-diagonal elements at column 0 and 2, values are -1 and -2)
1 1 -1           (Row 2: 1 off-diagonal element at column 1, value is -1)

Function to Implement

To help you build the CG algorithm, you need to implement the following helper functions:

• dot(v1, v2): Returns the dot product of two vectors v1 and v2(i.e.,,)
• vec_add(v1, v2): Returns a new vector which is the element-wise addition of v1 and v2.
• vec_sub(v1, v2): Returns a new vector which is the element-wise subtraction of v1 and v2.
• vec_const_mul(c, v): Returns a new vector where each element of v is multiplied by a constant c.
• spmul(d, pos, val, v): Preforms sparse matrix-vector multiplication and returns a new vector
• CG(A_d, A_pos, A_val, b, x0, step): The main Conjugate Gradient loop. Use the formulas provided above and your helper functions to find

Template code

def spmul(d:list,pos:list,val:list,v:list)->list:
    #TODO
    pass

def dot(v1:list,v2:list)->float:
    #TODO
    pass

def vec_add(v1:list,v2:list)->list:
    #TODO
    pass

def vec_sub(v1:list,v2:list)->list:
    #TODO
    pass

def vec_const_mul(c:float,v:list)->list:
    #TODO
    pass

def CG(A_d:list,A_pos:list,A_val:list,b:list,x0:list,step:int)->list:
    #TODO
    pass

import sys 
input = sys.stdin.readline 

if __name__=='__main__':
    n,step=map(int,input().split())
    b=list(map(int,input().split()))
    x0=list(map(int,input().split()))
    d=list(map(int,input().split()))
    pos=[ [] for _ in range(n) ]
    val=[ [] for _ in range(n) ]
    for i in range(n):
        tmp_list=list(map(int,input().split()))
        m=tmp_list[0]
        if m ==0 : continue
        tmp=tmp_list[1:]
        pos[i]=tmp[0:m]
        val[i]=tmp[m:2*m]
    
    ans_x = CG(d,pos,val,b,x0,step)
    for x in ans_x:
        print(f'{x:.5f}',end=' ')

Input

the first line contains two integers:
N step
The line contains integers, representing the vector b.
The line contains integers, representing the initial guess vector x0.
The line contains integers, representing d(the diagonal elements of ).
The next lines describe the off-diagonal non-zero elements for row to row . Each line follows this format:
M_i pos_i_0 pos_i_1 ... pos_i_{M-1} val_i_0 val_i_1 ... val_i_{M-1}

If the line will just be .

Constraints

• The martix will be a symmetric, positive-definite matrix.

Output

Output the final vector .
Print the elements of the vector separated by a single space. Please format each number to decimal places.

Sample Input  Download

Sample Output  Download

Discuss