15034 - Hard LA   

Description

You do not need to understand the following statement to solve the problem, you may skip directly to the Task section.

After repeatedly running out of paper, Cabi finally decided that relying on physical paper was fundamentally unsafe.

A sheet of paper may be lost, torn apart, given to suffering students during a lab, or destroyed by an unexpected earthquake.

Sometimes, Cabi falls down and causes a huge earthquake, which may destroy his remaining paper. Therefore, Cabi began designing an earthquake-resistant storage facility exclusively for his remaining paper.

Unfortunately, Cabi had no idea how to determine whether the structure would survive an earthquake.

To analyze the building, engineers model the displacement of a solid object using the equations of elastodynamics: , where

  • is the material density,
  • is the displacement field,
  • is the stress tensor,
  • represents external forces.

After studying for the equation for several minutes, Cabi successfully understood none of it ._. .

Fortunately, a mysterious engineer had already transformed the entire problem into a system of linear equations: .

The matrix is large, sparse, symmetric, and positive-definite.

Cabi first tried to solve the system using the Conjugate Gradient algorithm that you provided for him. However, the convergence was not fast enough to guarantee the safety of his paper before the next time he fell down.

He therefore decided to use a more powerful method: the Preconditioned Conjugate Gradient algorithm, also known as PCG.

Unfortunately, implementing it required paper to calculate the preconditioner.

Cabi had no paper.

Therefore, this problem now belongs to you.


TASK

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

Unlike Gaussian Elimination, CG is an iterative method perfect for large, sparse, and symmetric positive-definite matrices. However, there is an even stronger version of CG called Preconditioned Conjugate Gradient (PCG), which uses a preconditioner to improve the convergence of the algorithm.

Here is how it works:

  1. Prepare a starting vector called , the matrix and the vector .

  2. Initiate variables: let

    •  
    • is the diagonal of matrix , , and for all
    • , where
    •  
  3. Iterate the below for

    •  
    •  
    •  
    •  
    •  
    •  

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

Variables

  • N:The dimension of the linear system. This is also the size of the vectors and the matrix
  • step: the number of iterations you need to perform for the PCG 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 5 (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)

And for such it is

 

 

 

 

 

 

Function to Implement

To help you build the PCG 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): Performs sparse matrix-vector multiplication and returns a new vector
  • PCG(A_d, A_pos, A_val, b, x0, step): The main Preconditioned 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 PCG(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 = PCG(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 0.

Constraints

  •  
  •  
  •  
  •  
  •  
  •  
  •  
  • The matrix 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

Tags




Discuss