Givens求QR分解（3阶）_2017.2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2 2 1]\n",
      " [0 2 2]\n",
      " [2 1 2]]\n"
     ]
    }
   ],
   "source": [
    "#输入一个矩阵(3阶)\n",
    "A = numpy.array([[2, 2 , 1],[0, 2, 2],[2, 1, 2]])\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[2 0 2]\n"
     ]
    }
   ],
   "source": [
    "b1 = numpy.array([A[0,0],A[1,0],A[2,0]])\n",
    "print(b1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T121\n",
      "[[ 1.  0.  0.]\n",
      " [-0.  1.  0.]\n",
      " [ 0.  0.  1.]]\n",
      "T131\n",
      "[[ 0.70710678  0.          0.70710678]\n",
      " [ 0.          1.          0.        ]\n",
      " [-0.70710678  0.          0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "if A[0,0] != 0 or A[1,0] != 0 : \n",
    "    c1 = A[0,0]/(A[0,0]**2+A[1,0]**2)**0.5\n",
    "    s1 = A[1,0]/(A[0,0]**2+A[1,0]**2)**0.5\n",
    "    T121 = numpy.array([[c1, s1 , 0],[-s1, c1, 0],[0, 0, 1]])\n",
    "    print(\"T121\")\n",
    "    print(T121)\n",
    "    \n",
    "    c2 = (A[0,0]**2+A[1,0]**2)**0.5/(A[0,0]**2+A[1,0]**2+A[2,0]**2)**0.5\n",
    "    s2 = A[2,0]/(A[0,0]**2+A[1,0]**2+A[2,0]**2)**0.5\n",
    "    T131 = numpy.array([[c2, 0 , s2],[0, 1, 0],[-s2, 0, c2]])\n",
    "    print(\"T131\")\n",
    "    print(T131)\n",
    "    \n",
    "    \n",
    "elif (A[0,0]**2+A[1,0]**2)**0.5 != 0 or A[2,0] != 0 :\n",
    "    c2 = (A[0,0]**2+A[1,0]**2)**0.5/(A[0,0]**2+A[1,0]**2+A[2,0]**2)**0.5\n",
    "    s2 = A[2,0]/(A[0,0]**2+A[1,0]**2+A[2,0]**2)**0.5\n",
    "    T131 = numpy.array([[c2, 0 , s2],[0, 1, 0],[-s2, 0, c2]])\n",
    "    print(\"T131\")\n",
    "    print(T131)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T121 b1\n",
      "[2. 0. 2.]\n",
      "T131(T121 b1)\n",
      "[2.82842712 0.         0.        ]\n",
      "T1\n",
      "[[ 0.70710678  0.          0.70710678]\n",
      " [ 0.          1.          0.        ]\n",
      " [-0.70710678  0.          0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "if A[0,0] != 0 or A[1,0] != 0 : \n",
    "    \n",
    "    print(\"T121 b1\")\n",
    "    print(numpy.dot(T121 ,b1))\n",
    "    print(\"T131(T121 b1)\")\n",
    "    print(numpy.dot(T131,numpy.dot(T121,b1)))    \n",
    "    T1 = numpy.dot(T121,T131)      \n",
    "    print(\"T1\")    \n",
    "    print(T1)\n",
    "    \n",
    "    \n",
    "elif(A[0,0]**2+A[1,0]**2)**0.5 != 0 or A[2,0] != 0 :\n",
    "   \n",
    "    T1 = T131\n",
    "    print(\"T1\")    \n",
    "    print(T1)\n",
    "    print(\"T131(T121 b1)\")\n",
    "    print(numpy.dot(T1,b1))\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T1 A\n",
      "[[ 2.82842712  2.12132034  2.12132034]\n",
      " [ 0.          2.          2.        ]\n",
      " [ 0.         -0.70710678  0.70710678]]\n",
      "A1\n",
      "[[ 2.          2.        ]\n",
      " [-0.70710678  0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "M = numpy.dot(T1 ,A)\n",
    "print(\"T1 A\")\n",
    "print(M)\n",
    "print(\"A1\")\n",
    "A1 = numpy.array([[M[1,1],M[1,2]],[M[2,1],M[2,2]]])\n",
    "print(A1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 2.        ]\n",
      " [-0.70710678]]\n"
     ]
    }
   ],
   "source": [
    "b2 = numpy.array([[A1[0,0]],[A1[1,0]]])\n",
    "print(b2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T122\n",
      "[[ 0.94280904 -0.33333333]\n",
      " [ 0.33333333  0.94280904]]\n",
      "T122b2\n",
      "[[2.12132034]\n",
      " [0.        ]]\n",
      "T2\n",
      "[[ 0.94280904 -0.33333333]\n",
      " [ 0.33333333  0.94280904]]\n"
     ]
    }
   ],
   "source": [
    "if A1[0,0] != 0 or A1[1,0] != 0 : \n",
    "    c1 = A1[0,0]/(A1[0,0]**2+A1[1,0]**2)**0.5\n",
    "    s1 = A1[1,0]/(A1[0,0]**2+A1[1,0]**2)**0.5\n",
    "    T122 = numpy.array([[c1, s1],[-s1, c1]])\n",
    "    print(\"T122\")\n",
    "    print(T122)\n",
    "\n",
    "    T122b2 = numpy.dot(T122 ,b2)\n",
    "    print(\"T122b2\")\n",
    "    print(T122b2)\n",
    "    \n",
    "    T2 = T122\n",
    "    print(\"T2\")\n",
    "    print(T2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "T2 A1\n",
      "[[2.12132034 1.64991582]\n",
      " [0.         1.33333333]]\n"
     ]
    }
   ],
   "source": [
    "N = numpy.dot(T2 ,A1)\n",
    "print(\"T2 A1\")\n",
    "print(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.70710678  0.          0.70710678]\n",
      " [ 0.23570226  0.94280904 -0.23570226]\n",
      " [-0.66666667  0.33333333  0.66666667]]\n"
     ]
    }
   ],
   "source": [
    "T = numpy.dot([[1, 0, 0],[0,T2[0,0],T2[0,1]],[0,T2[1,0],T2[1,1]]],T1)\n",
    "print(T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.70710678  0.23570226 -0.66666667]\n",
      " [ 0.          0.94280904  0.33333333]\n",
      " [ 0.70710678 -0.23570226  0.66666667]]\n"
     ]
    }
   ],
   "source": [
    "Q = T.T\n",
    "print(Q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2.82842712 2.12132034 2.12132034]\n",
      " [0.         2.12132034 1.64991582]\n",
      " [0.         0.         1.33333333]]\n"
     ]
    }
   ],
   "source": [
    "R = numpy.array([[M[0,0],M[0,1],M[0,2]],[0, N[0,0],N[0,1]],[0, N[1,0],N[1,1]]])\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[2. 2. 1.]\n",
      " [0. 2. 2.]\n",
      " [2. 1. 2.]]\n"
     ]
    }
   ],
   "source": [
    "A = numpy.dot(Q,R)\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}