1// Copyright ©2015 The Gonum Authors. All rights reserved. 2// Use of this source code is governed by a BSD-style 3// license that can be found in the LICENSE file. 4// 5// Some of the loop unrolling code is copied from: 6// http://golang.org/src/math/big/arith_amd64.s 7// which is distributed under these terms: 8// 9// Copyright (c) 2012 The Go Authors. All rights reserved. 10// 11// Redistribution and use in source and binary forms, with or without 12// modification, are permitted provided that the following conditions are 13// met: 14// 15// * Redistributions of source code must retain the above copyright 16// notice, this list of conditions and the following disclaimer. 17// * Redistributions in binary form must reproduce the above 18// copyright notice, this list of conditions and the following disclaimer 19// in the documentation and/or other materials provided with the 20// distribution. 21// * Neither the name of Google Inc. nor the names of its 22// contributors may be used to endorse or promote products derived from 23// this software without specific prior written permission. 24// 25// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 26// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 27// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 28// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 29// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 30// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 31// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 32// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 33// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 34// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 35// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 36 37// +build !noasm,!appengine,!safe 38 39#include "textflag.h" 40 41#define X_PTR SI 42#define Y_PTR DI 43#define DST_PTR DI 44#define IDX AX 45#define LEN CX 46#define TAIL BX 47#define ALPHA X0 48#define ALPHA_2 X1 49 50// func AxpyUnitary(alpha float64, x, y []float64) 51TEXT ·AxpyUnitary(SB), NOSPLIT, $0 52 MOVQ x_base+8(FP), X_PTR // X_PTR := &x 53 MOVQ y_base+32(FP), Y_PTR // Y_PTR := &y 54 MOVQ x_len+16(FP), LEN // LEN = min( len(x), len(y) ) 55 CMPQ y_len+40(FP), LEN 56 CMOVQLE y_len+40(FP), LEN 57 CMPQ LEN, $0 // if LEN == 0 { return } 58 JE end 59 XORQ IDX, IDX 60 MOVSD alpha+0(FP), ALPHA // ALPHA := { alpha, alpha } 61 SHUFPD $0, ALPHA, ALPHA 62 MOVUPS ALPHA, ALPHA_2 // ALPHA_2 := ALPHA for pipelining 63 MOVQ Y_PTR, TAIL // Check memory alignment 64 ANDQ $15, TAIL // TAIL = &y % 16 65 JZ no_trim // if TAIL == 0 { goto no_trim } 66 67 // Align on 16-byte boundary 68 MOVSD (X_PTR), X2 // X2 := x[0] 69 MULSD ALPHA, X2 // X2 *= a 70 ADDSD (Y_PTR), X2 // X2 += y[0] 71 MOVSD X2, (DST_PTR) // y[0] = X2 72 INCQ IDX // i++ 73 DECQ LEN // LEN-- 74 JZ end // if LEN == 0 { return } 75 76no_trim: 77 MOVQ LEN, TAIL 78 ANDQ $7, TAIL // TAIL := n % 8 79 SHRQ $3, LEN // LEN = floor( n / 8 ) 80 JZ tail_start // if LEN == 0 { goto tail2_start } 81 82loop: // do { 83 // y[i] += alpha * x[i] unrolled 8x. 84 MOVUPS (X_PTR)(IDX*8), X2 // X_i = x[i] 85 MOVUPS 16(X_PTR)(IDX*8), X3 86 MOVUPS 32(X_PTR)(IDX*8), X4 87 MOVUPS 48(X_PTR)(IDX*8), X5 88 89 MULPD ALPHA, X2 // X_i *= a 90 MULPD ALPHA_2, X3 91 MULPD ALPHA, X4 92 MULPD ALPHA_2, X5 93 94 ADDPD (Y_PTR)(IDX*8), X2 // X_i += y[i] 95 ADDPD 16(Y_PTR)(IDX*8), X3 96 ADDPD 32(Y_PTR)(IDX*8), X4 97 ADDPD 48(Y_PTR)(IDX*8), X5 98 99 MOVUPS X2, (DST_PTR)(IDX*8) // y[i] = X_i 100 MOVUPS X3, 16(DST_PTR)(IDX*8) 101 MOVUPS X4, 32(DST_PTR)(IDX*8) 102 MOVUPS X5, 48(DST_PTR)(IDX*8) 103 104 ADDQ $8, IDX // i += 8 105 DECQ LEN 106 JNZ loop // } while --LEN > 0 107 CMPQ TAIL, $0 // if TAIL == 0 { return } 108 JE end 109 110tail_start: // Reset loop registers 111 MOVQ TAIL, LEN // Loop counter: LEN = TAIL 112 SHRQ $1, LEN // LEN = floor( TAIL / 2 ) 113 JZ tail_one // if TAIL == 0 { goto tail } 114 115tail_two: // do { 116 MOVUPS (X_PTR)(IDX*8), X2 // X2 = x[i] 117 MULPD ALPHA, X2 // X2 *= a 118 ADDPD (Y_PTR)(IDX*8), X2 // X2 += y[i] 119 MOVUPS X2, (DST_PTR)(IDX*8) // y[i] = X2 120 ADDQ $2, IDX // i += 2 121 DECQ LEN 122 JNZ tail_two // } while --LEN > 0 123 124 ANDQ $1, TAIL 125 JZ end // if TAIL == 0 { goto end } 126 127tail_one: 128 MOVSD (X_PTR)(IDX*8), X2 // X2 = x[i] 129 MULSD ALPHA, X2 // X2 *= a 130 ADDSD (Y_PTR)(IDX*8), X2 // X2 += y[i] 131 MOVSD X2, (DST_PTR)(IDX*8) // y[i] = X2 132 133end: 134 RET 135